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Abstract 


The similarity principle is employed to determine the maximum velocity of the vertical 
propagation of a diffusing smoke flowing from a continuous fixed source in a stratified fluid. 
With the use of the velocity profile (Kao, 1959) for a stratified fluid the profile of the boundaries 
of the smoke plume is evaluated numerically for various values of the stability parameter 


“= kgzuQ (c, QT) ‘ui. It is shown that the spreading of a smoke plume generally in- 


creases with increasing heat flux from the boundary surface. 


I. Introduction 


Since the publication of Taylor’s pioneering 
paper, “Diffusion by continuous movements” 
(TayLoR, 1921), in 1921, the statistical charac- 
teristics of turbulent diffusion of particles 
released serially from a fixed point in a steady, 
homogenous and stratified fluid have been 
studied by many investigators, notably by 
RICHARDSON (1926), SUTTON (1932), BATCHE- 
LOR (1949, 1950, 1952), BRIER (1950), EDINGER 
(1953), FRENKIEL (1953), Panorsky and Mc- 
Cormick (1954), CORRSIN (1956), ELLISON 
(1957), GIFFORD (1957, 1959), TOWNSEND 
(1958), Barap (1959), INOUE (1959), and 
STEWART (1959). Advances in both experi- 
ments and theories have been made. 

In an analysis of the effects of stratification 
on the turbulent transfers in the steady atmos- 
pheric flow over an infinite uniform plane 
Kao (1959) has shown that for a constant 
free stream velocity the mean velocity in the 
surface boundary layer may be expressed in 
terms of the dimensionless height, 7 =z/z, and 


1 Present affiliation: Department of Meteorology, 
University of Utah, Salt Lake City, Utah. 
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the dimensionless parameters sy) and f defined 


respectively by 
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net (x) 
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where k is von Karman constant, g is the 
gravity acceleration, Q is the vertical heat 
flux, 2, is the roughness height, c, is the specific 
heat of the air at constant pressure, u,, is the 


friction velocity in neutral stratification, B is 
a constant, H, is the height of the free stream 


velocity, @ and T are respectively the mean 
density and temperature of the air. Applying 
the similarity principle to the study of the 
turbulent flow in the surface layer of a stratified 
fluid Monin (1959) has analysed the character- 
istics of turbulent diffusion and obtained some 
empirical results. The objective of this paper 
is to give a theoretical estimate of the maximum 
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velocity of the vertical propagation of the 
diffusing smoke and the profile of the bounda- 
ries of the smoke plume making use of the 
above mentioned velocity profile for a strati- 


fied fluid. 


2. Energy Balance Equation and the Maxi- 
mum Velocity of the Vertical Propaga- 
tion of the Diffusing Smoke 


The Navier-Stokes equations of motion and 
the equation of continuity for a viscous in- 
compressible fluid may respectively be written 


EA eee ae 
Con CEE CEE 
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where repeated suffixes imply summation over 
three values of the repeated suffix, p is the 
pressure, u; and — g; are respectively the i-com- 
ponents of the velocity and the acceleration of 
the gravitational field. The barred and primed 
quantities in (3) and (4) represent respectively 
the time-averaged and turbulent quantities 


defined as 


t+T 
I 
ana) f obnzdd (5) 
t-T 
pg=p-p (6) 


Multiplying (3) by w; then taking time- 
average of each member of the equation, we 
obtain the equation for the turbulent kinetic 
energy per unit mass 
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where - g?= - (u? + ui? + ur?) is the kinetic 
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energy of the turbulent flow. On the left 
hand side of (7) the first term is the local 
rate of change of kinetic energy of the tur- 
bulent flow, the second term represents the 
transport of kinetic energy by diffuse motion, 
the third term represents the transport of 
kinetic energy through advection by the mean 
motion. The remaining terms in (7) represent 
energy production by transfer from the mean 
flow and energy loss through work done 
against buoyancy forces and through viscous 
dissipation. 

Let the x-axis be directed along the mean 
velocity and the z-axis vertical. Taking average 


of (7) over the xy-plane, defined by 


7,1; 7)= 
x Y 
= lim rl G(x, y,z,t;t)dxdy (8) 
Rey 
we obtain 
een 
Ot 21 oz 


(9) 


Here equation of continuity is employed, 


u = u and T =T are assumed. 

It is well-known that for high Reynolds 
numbers the rates of dissipation are very 
largely determined by the typical length and 
velocity of the turbulence and not by molec- 
ular quantities such as viscosity. In accordance 
with the similarity principle the rate of dissipa- 
tion of turbulent energy may be expressed by 
(w’)8/2]=1 where | is the scale of turbulence. 
Introducing the concepts of eddy viscosity Kj, 
and eddy conductivity Kj, defined respectively 
by the following momentum- and heat-transfer 
equation 


du = 
1e RP UT uw 
and 
do =, 
KT =-wT 
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in (9), we obtain for the stationary case the 
energy balance equation 


du\? 
= Rn (&) Fam K;, 


where wy -(w?)}. Since K„aw,l, (10) may 
be written 


dO £ 
ET (10) 


er: (I - R,)* (11) 
th 
where 
g4 
Ky "T4 
IF Tg (12) 


is the flux Richardson number. 

According to the similarity principle the 
maximum velocity of the vertical propagation 
of the diffusing smoke may be written 


(13) 


where I(syn) is a universal function, and is 
subjected to the condition that I> 1 as |s9| +0. 
B is therefore equal to wy/u, under neutral 
stratification. MoNIN (1959) obtained from ex- 
periment that $ ~ 0.75. 

It has been shown (Kao, 1959) that 
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Substitution of (13) and (14) in (11) gives 
S 4 
pe 

; E + (: + ss n)'| 


i ig 


For small values of son, (15) may be approxi- 
mated by 


(15) 


I(son)=!ı+ 


I 


«y 
I(son)¥1+—- 59 (16) 
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The ratio of the maximum velocity of the 
vertical propagation of the smoke and the 
friction velocity increases with son. 


3. The Shape of the Boundaries of the Smoke 
Plume 


We consider a smoke plume of neutral 
stratification flowing out of a stationary point 
source at the height z, in the surface layer 
of the stratified fluid. The equations for smoke 
particles at the boundaries of the plume may 
be written 


(17) 


dm = 


where the positive and negative signs refer 
respectively to the upper and lower boundaries 
of the smoke plume. In view of (13) and (15) 
the equations for the boundaries of the smoke 
plume take the form 


dx 


G(on) 
De ae 


I(son) 


where G(so7) has been shown (Kao, 1959) 
to be 


(18) 
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Table 1. Values of & (son) 


100 800 900 1000 


et 
5 |14.007|358.499)847-386|1380.547|1939.017[2513-432|3097-852|3687.705|4278.283|4865-5325458-419 
4 14.009]358.861|848.991|1384-404|1946.263|2525-368]3116.036/3714-174|4316.344/4918.823)5517.163 
3 1.4.012|359.319|851.017|1389.256|1955.327|2540-187|3138.362|3746.075|4360.444|/4979-043|5599-531 
2 
I 


300 400 500 


600 | 700 


14.016/359.895|853.556|1395-310|1966.577|2558-450[3165.613|3784-492|4412.480]5047.558|5688.080 
14.020|360.614|856.721|1402.823|1980.465|2580.850|3198.775|3830.787'4474-404|5127:776|5789-443 

o 14.026|361.517|860.663|1412.134|1997.584|2608.300|3236.153/3886.750/4548.68015223.143|5908.740 
— I 1 4.033|362-642|865.556/1423-639|2018.645|/2641.925|3288.395|3954-701|4038.461|5337-899|6051.637 
— 2 1 4.042]364.039|871.620|1437-847|2044.574|2683.212|3348.723|4037-796|4748.088]5477-865|/6225-792 
== 5) 14.052|365.784|870.153|1455.452/2076.643|2734.234|3423.282|4140.581|4883.905|5651.650|6442.629 
—4 14.067|367.966|888.539|1477-353\2116.559|2797-868|3516.578|4269.768|5055.547|5872-703|6720.501 
—5 14.085/370.704|900.296|1504.792|2166.755|2878.364|3035-524|4436.006|5279.036|6164.377|7092-789 


For small values of syn, the diffusion angle 7 


b imated b HN = pee G (sw) ge 
may be approximated by E (son; ns) = PB = & Ts) 
Pets a 
"3 
an € ~ tan 1} Bk . (20) Here de 
ungen me 
Integration of (18) from the dimensionless “G(s t) 
source-height, 7,=z,/Z to an arbitrary height = an de} = 
n gives the profile of the boundaries of the : (Soû 
smoke plume as a function of dimensionless 
height n and stability parameter 5, = + {E (son) — E(sons)} 
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Fig. 1. The variation of £ as a function of n, for various values of the 
stability parameter sp. 
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where 


Integral (22) has been evaluated numerically 
for various values of sy and 7, and are shown 
in Table 1 in which B=1 and H,/z,=3 - 10° are 
chosen. It is seen from Table 1 that for a 
fixed value of 7 the value of & (sox) decreases 
with increasing value of sy. This indicates that 
the spreading of a smoke plume increases with 
increasing value of sy. It may be noted that 
at 7 = 10° the difference in the values of & for 
s=-5-ro”"* and the neutral case is more 
than 17% of the value for the latter. The 
E (son)-curves for various values of sy and 7 
are shown in Fig. 1. 


4. Conclusion 
The effects of thermal stratification on tur- 
bulent diffusion from a continuous fixed source 
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are examined. The similarity principle is em- 
ployed to determine the maximum velocity of 
propagation of a diffusing smoke in a stratified 
fluid. The profile of the boundaries of the 
smoke plume is evaluated numerically for 
various values of the stability parameter 
So = kgzoQ(c OT) ug}. It is shown that the 
spreading of a smoke plume generally in- 
creases with increasing heat flux from the 
boundary surface. 
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Energy and Numerical Weather Prediction 
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Abstract 


Since the study of energy transformations and the numerical integration of simplifiied equa- 
tions are sometimes used as alternative approaches to the same physical problem, it is often 
desirable that the simplified equations conserve total energy under reversible adiabatic pro- 
cesses. Preferably, the equations should also conserve the sum of kinetic energy and available 
potential energy, and they should describe the tendency for static stability to increase as kinetic 
energy is released. 

It is found that if the equation of balance is used as a filtering approximation, all the terms in 
the vorticity equation which involve both the rotational and the divergent part of the wind 
field should be retained, while, if the geostrophic equation is used, all of these terms in the 
vorticity equation should be omitted, if the equations are to possess suitable energy invariants. 

An n-layer model with the appropriate energy invariants is developed. The two-layer model 
may be the simplest possible model with variable static stability. The model appears to be 
suitable for theoretical studies of the general circulation. 


I. Introduction 


One enlightening method of studying the 
behavior of the atmosphere, or a portion of it, 
consists of examining the behavior of the 
energy involved. Any atmospheric circulation 
system, whether it be a small-scale convection 
cell, a cyclone, or a large-scale zonal-wind 
system, is marked by a supply of kinetic 
energy, and the development of such a system 
requires either a transformation of some 
other form of energy into kinetic energy, or a 
conversion of the kinetic energy of some other 
system into that of the developing system. 
The classical paper adopting this method is 
that of MARGULES (1903); more recently nu- 
merous papers concerning energy transforma- 
tions in the atmosphere have greatly increased 


! The research reported in this work has been sponsored 
by the Geophysics Research Directorate of the Air 
Force Cambridge Research Center, under Contract No. 
AF 19(604)-1000. 


our understanding of the general circula- 
tion. 

Another method which is currently finding 
much favor consists of generating sequences of 
“weather maps”, through numerical integra- 
tion of a simplified set of dynamic equations— 
ordinarily a set which has been developed for 
use in numerical weather prediction. A paper 
of Pmurrips (1956), which aptly demonstrates 
the power of this method, has already become 
a classic. 

There are a number of problems to which 
both of these methods are applicable. If the 
methods are to yield compatible results, it 
would appear desirable that the simplified 
equations should not violate any of the energy 
principles which the exact equations express. 
Hence even the simplified equations should 
conserve total energy under reversible adiabatic 
processes, and they should lead to approriate 
expressions for the conversion of one form of 
energy into another. 
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2. Energy, available potential energy, and 
gross static stability 


Of the various forms of energy present in 
the atmosphere, kinetic energy has often 
received the most attention. Often the total 
kinetic energy of a weather system is regarded 
as a measure of its intensity. The only other 
forms of atmospheric energy which appear 
to play a major role in the kinetic energy 
budget of the troposhere and lower stratosphere 
are potential energy, internal energy, and the 
latent energy of water vapor. Potential and 
internal energy may be transformed directly 
into kinetic energy, while latent energy may 
be transformed directly into internal energy, 
which is then transformed into kinetic energy. 

It is easily shown by means of the hydro- 
static approximation that the changes of the 
potential energy P and the internal energy J of 
the whole atmosphere are approximately pro- 
portional, so that it is convenient to regard 
potential and internal energy as constituting 
a single form of energy. This form has been 
called total potential energy by MARGULES (1903). 

In terms of potential temperature ©, the 
total potential energy of the whole atmosphere 
may be given by 


P+I=c,po "JS p* OGM, (1) 


where p = pressure, pop is a standard pressure 
(1,000 mb), x is the ratio (c,—c,)/c,, c, and 
cp are the specific heats of air, and dM is an 
element of mass of the atmosphere. 

In the long run, there must be a net depletion 
of kinetic energy by dissipative processes. It 
follows that there must be an equal net 
generation of kinetic energy by reversible 
adiabatic processes; this generation must occur 
at the expense of total potential energy. It 
follows in turn that there must be an equal net 
generation of total potential energy by heating 
of all kinds. These three steps comprise the 
basic energy cycle of the atmosphere. The 
rate at which these steps proceed is a fundamen- 
tal characteristic of the general circulation. 

The writer (1955, 1960) has shown that a 
partial explanation of the intensity of the 
energy cycle can be obtained by considering 
available potential energy. Available potential 
energy does not represent a supply of energy 
additional to the forms already mentioned, 
but instead represents a portion of the total 
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potential energy which may be available for 
conversion into kinetic energy. It is equal to 
the excess of total potential energy, above the 
total potential energy which would be present 
if the mass of the atmosphere were to be re- 
arranged, under isentropic changes of state, 
to possess horizontal isentropic surfaces, with 
stable stratification. 

Under this hypothetical rearrangement of 
mass, @ and p completely determine cach 
other. Since this rearrangement would conserve 
the total mass lying above a given isentropic 
surface, it would conserve the average value p 
of p on each isentropic surface, and the resulting 
value of P+ I would be obtained by replacing 
p by p in (4). Hence the available potential 
energy is given by 


A=cpse f (p* — Pr) OdM. (2) 


The writer (1955) has derived from (2) the 
approximate expression 


I Sef Ys, TONER vs 
A~ Exopie |p (=) (© - @)?dM, (3) 


where a bar (~) denotes an average over an 
isobaric surface. Thus A is approximated by a 
weighted average of the variance of © within 
isobaric surfaces. 

It would appear desirable, then, that in any 
numerical study of the energetics of the 
atmosphere, the equations used should conserve 
the sum of kinetic energy and available 
potential energy, under reversible adiabatic 
processes. 

Another quantity which plays a part in the 
energetics of the atmosphere is static stability. 
As its name might imply, static stability has 
long been regarded as indicating the tendency 
for convective overturning to develop in an 
atmosphere at rest. More recently, it has been 
recognized as a factor in determining the 
dynamic stability of a baroclinic flow. 

The generation of kinetic energy appears 
to be accompanied by a sinking of colder air 
and a simultaneous rising of warmer air across 
the same levels. (It must be accompanied by a 
pressure increase of the colder air and a pressure 
decrease of the warmer air.) This process 
should, by lifting the warmer air above the 
colder air, increase the static stability in some 
over-all sense. 
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In a manner analogous to the definition of 
available potential energy, we can define a 
quantity, which we shall call gross static stability, 
whose variations under reversible adiabatic 
processes are equal to those of kinetic energy. 
Gross static stability is equal to the deficit of 
total potential energy, below the total poten- 
tial energy which would be present if the 
mass of the atmosphere were to be rearranged, 
under isentropic changes of state, to possess 
vertical isentropic surfaces. 

Under this hypothetical rearrangement of 
mass, © and p are completely independent. The 
resulting value of P + J after this rearrange- 
ment would therefore be obtained by re- 
placing p* in (1) by its average value over the 
mass of a vertical column, 1.e., by (I + x)-1p# 
where po is the surface pressure. Hence the 
gross static stability is given by 


so that S is expressible as a weighted average 
of —0O0/dp, which may be taken as a measure 
of static stability. 

It would thus also appear desirable that in a 
numerical study of the energetics of the at- 
mosphere, the equations used should conserve 
the difference between kinetic energy and 
gross static stability, under reversible adiabatic 
processes. Equations which allow the release 
of kinetic energy without static stabilization 
could conceivably overpredict the growth of 
disturbances. 


3. Energy and the filtering approximations 


In order to assess the suitability of the equa- 
tions of numerical weather prediction for 
studying the energetics of the atmosphere, we 
shall first outline a procedure for obtaining 
these equations from the equations which 
directly express the laws governing the atmos- 
sphere. It might be noted that historically 
the development of the simplified equations 
has reached its present position by a somewhat 
different route. 


EDWARD N. LORENZ 


For a dry atmosphere, the physical laws 
determine a set of five scalar prognostic equa- 
tions—three equations of motion, the equation 
of continuity, and the thermal equation—and 
one nonprognostic equation or identity—the 
equation of state. These equations contain six 
dependent variables; the prognostic equations 
may be expressed in terms of five dependent 
variables with the aid of the one identity. 

The equation of vertical motion is first 
discarded, and replaced by the hydrostatic 
equation—and identity. The equation of con- 
tinuity and the thermal equation reduce to one 
prognostic equation and one identity with the 
aid of the time derivative of the hydrostatic 
equation. Thus there remain three prognostic 
equations, which may be expressed in terms 
of three dependent variables with the aid of 
the three identities. This system of equations 
is the one generally called the primitive equa- 
tions. 

The new system is next expressed with 
pressure as an independent variable, and height 
as a dependent variable. The horizontal wind 
components are then expressed in terms of 
vorticity and divergence, and the equations 
of horizontal motion are expressed by their 
equivalents—the vorticity equation and the 
divergence equation. 

The divergence equation is then discarded, 
and replaced by the equation of balance, an 
identity, obtained by dropping from the diver- 
gence equation all the terms which contain 
divergence. The vorticity equation and the 
thermal equation reduce to one prognostic 
equation and one identity with the aid of the 
time derivative of the equation of balance. 
Thus there remains one prognostic equation, 
which may be expressed in terms of one de- 
pendent variable with the aid of the five 
identities. 

It is often more convenient to omit certain 
additional terms from the equation of balance, 
reducing it to the geostrophic equation. Cer- 
tain terms in the vorticity equation are also 
often omitted. The new system still contains 
but one prognostic equation. 

The equation of balance and the geostrophic 
equation are often called filtering approxima- 
tions, since they eliminate the occurrence of 
certain waves which are often considered 
irrelevant and which can occur in systems 
governed by the primitive equations. 
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Finally, the vertical dimension may be 
replaced by several layers. Each function of 
time and three space dimensions is then re- 
placed by several functions of time and two 
space dimensions. 

With the original set of five prognostic 
equations as the governing equations, total 
energy is conserved under reversible adiabatic 
processes. After the equation of vertical motion 
is replaced by the hydrostatic equation, total 
energy may still be said to be conserved, but 
only if the kinetic energy contained in the 
vertical component of the motion is not in- 
cluded in the total amount of kinetic energy. 
Since the omitted kinetic energy is an insig- 
nificant fraction of the total (provided that 
we are dealing with large-scale systems), this 
restriction is of little consequence. 

Let us see what happens when further modi- 
fications are made. Choosing pressure as the 
vertical coordinate, let the horizontal wind V 
be written as 


Vi Ve Ve (6) 


where V, is nondivergent and V3 is irrotational. 
We shall attach a subscript “2” or “3” to any 
dependent variable related to V, or V; through 
an identity. Thus we may introduce a stream 
function y, and a velocity potential y, such 
that 


Va: =k X VWo (7) 
Vs = VX (8) 


where k is the vertical unit vector. 
The vorticity €, and the divergence 6, then 
satisfy the relations 


C.=v-Vxk=v-V,xk=v’y,, (9) 
(10) 


while the individual pressure change w 3 is 
related to 6, through the equation of conti- 


nuity G ra 


Likewise we shall attach a subscript “1” to 
the temperature T, and to any dependent 
variable related to it through an identity. Thus 
specific volume &, and elevation z, appear in 
the equation of state 


63 =V-V=v-V3=V9743 


d, + Jo,/p =0 


(12) 


p&ı = Xl 
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and the hydrostatic equation 


d21/0p + «lg =0, (13) 
where g is the acceleration of gravity, while 
potential temperature, @,, is defined by the 
relation 

(14) 


Equations (7) through (14) enable us to 
express any variable with a numerical sub- 
script as a linear function of any other variable 
with the same subscript. Hence any three 
variables with three different subscripts may 
be regarded as the dependent variables in the 
prognostic equations. 

The three prognostic equations, namely the 
thermal equation, the vorticity equation, and 
the divergence equation, may now be written 


Fp = Sa: 9) - Ver, (15) 


9, Fon el 


ae 
an = - (2S) —J (ve f) = ¥-f V3 - 
= Mae Woo = La 02 marge — V@3° 
pr _ I%s 
dp (os ine) 
003 


oV 
- v-(V;-vV;) VO3 * ap 
oN 
- v-(V3-4V5) re (17) 


provided that nonadiabatic effects are omitted. 
Here J denotes a Jacobian, e.g., 


J (%e; 0,) = VYa- VO, x k 


and f is the Coriolis parameter. 

The proper lower boundary conditions are 
z =o and dz/dt =o, if the earth’s surface is 
assumed horizontal. It is convenient at this 
point to simplify the system of equations by 
discarding these boundary conditions, and 
replacing them by the conditions p = p 9 = con- 
stant and w=o. With these new boundary 
conditions, total energy will still be conserved, 
while, since the statistical distribution of © will 


(18) 
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be conserved, the relations involving available 

otential energy and gross static stability will 
still hold. The new boundary conditions do 
not assume a flat sea-level pressure field, since 
z is no longer assumed constant at the lower 
boundary. 

From identities (6), (7), and (8), it follows 
that the kinetic energy per unit mass is 


it I 
VV => ye Ve +I (War Xe) + 


(19) 


I 
ye ie VYs- 


Since the average value of a Jacobian through- 
out the atmosphere vanishes, the kinetic 
energy of the whole atmosphere is given by 


We have already seen that the total potential 
energy of the whole atmosphere may be given 


by 
(21) 


The terms in the divergence equation (17) 
may be grouped into six classes, such that the 
different terms in any one class contain the 
same set of numerical subscripts. Thus the six 
classes may be denoted by (1), (2), (3), (2,2), 
(2,3) and (3,3). Likewise, the terms of the 
vorticity equation (16) fall into the five 
classes (2), (3), (2,2), (2,3) and (3,3), while 
those on the right of the thermal equation (15) 
fall into the two classes (1, 2) and (1,3). 

From the prognostic equations we may 
determine the classes into which the various 
terms fall, in the expressions for 2(P, +1,)/ot, 
OK,/dt, and 9K;/dt. In determining these 
classes, we shall make repeated use of integra- 
tion by parts, and observe that the divergence 
of any vector, the Jacobian of any two scalars, 
and the vertical derivative of any quantity 
which vanishes at the bottom and the top of 
the atmosphere, all vanish when averaged 
throughout the atmosphere. 

We then find that the only nonvanishing 
terms of 9(P,+1,)/dt fall into the class (1,3) 


P, Ar Ik, = CoP eat J p*O,dM. 
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while the nonvanishing terms of 2K,/0t fall 
into the classes (2,3), (2,2,3), and (2,3,3), and 
those of 9K;/dt fall into the classes (1,3), (2,3), 
(2,2,3), and (2,3,3). In particular, 


2 (Pi +1)=xi,p5* fp" 1O:wsdp. (22) 
Since equations (15), (16), and (17) conserve 
total energy, the terms of 2K;/9f in class (1,3) 
must cancel the expression for 2(P, +1,)/ot. 
The remaining terms of JK;/dt must then 
cancel, class by class, the terms of JK,/dt. 

Now consider what happens when the 
system. of equations is simplified by omitting 
certain terms from the divergence equation 
(17). First, if the term 90,/0f is omitted, (17) 
becomes an identity, and the statement that 
O(P, +1, + K,+Ks)dt vanishes must be re- 
placed by the statement that 9(P, +, +Kg) /ot 
vanishes. Thus the new system may still be 
said to preserve total energy, but only if Kg 
is not included in the total amount of kinetic 
energy. This restriction is quite analogous to 
the exclusion of the kinetic energy contained 
in the vertical motion, when the hydrostatic 
equation is first introduced. 

If the remaining terms in (17) which contain 
a subscript “3” are omitted, (17) reduces to 
the equation of balance. This omission results 
in the omission of the terms of class (2,3,3) 
from the expression for 7K,/dt. In order that 
total energy be still conserved, the terms of 
class (2,3,3) must be omitted from the equa- 
tion for JK,/dt, which is accomplished by 
omitting the terms of class (3,3) from the 
vorticity equation (16). In most previous 
studies these terms have been omitted as a 
matter of course. 

Further simplifications result from omitting 
the terms of class (2,2) from the divergence 
equation (17) (which has already been reduced 
to the equation of balance). The equation then 
becomes a form of the geostrophic equation 
This omission results in the omission of terms 
of class (2,2,3) from 9K,/0f. In order that the 
new system of equations may preserve total 
energy, it is thus necessary to omit the terms 
of class (2,2,3) from the equation for K,/dt, 
which is accomplished by omitting the terms 
of class (2,3) from the vorticity equation (16). 

In previous studies the four terms of class 
(2,3) in (16) have often, but by no means 
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invariably, been omitted. Of these four terms 
the second and third, which represent the 
change of vorticity due to concentration of 
contrasting currents, and the vertical advection 
of vorticity, have most frequently been in- 
cluded. The fourth term, often called the 
twisting term, may be equally important, and, 
as shown by REED and SANDERS (17), may be 
included with little additional difficulty. The 
first of these four terms, the advection of 
vorticity by the divergent part of the wind, 
seems to have been generally neglected. To the 
writer this neglect seems somewhat illogical 
when the other three terms are included; the 
presence of any vertical flow, which may 
advect vorticity, implies by continuity the 
presence of divergent horizontal flow, which 
may also advect vorticity. 

It now appears that all four of these terms 
should be included if the equation of balance 
is to be used, and all should be omitted if the 
geostrophic equation is to be used, in any 
study where the energetics are important. 
The inclusion of these terms, together with 
the geostrophic equation, or the omission of 
these terms, together with the equation of 
balance, yields a system of equations without a 
suitable energy invariant. 

According to equations (2) and (3), the 
sum of P+J and S will be conserved by any 
equations which conserve the average value 
of ©, while the difference of P+ I and A will 
also be conserved if the entire statistical distri- 
bution of © is conserved. The difference of S 
and K, and the sum of A and K, will then be 
conserved if the sum of P+I and K is con- 
served. Since we have not yet tampered with 
the thermal equation (15), these conditions 
are still satisfied. 

We must note that the term -V,:VO, 
representing advection by the divergent part 
of the wind, has not been omitted from (15). 
Like the advection of vorticity by V3, this 
term has been neglected in many studies. If 
the only modification of (15) is the omission of 
this term, the equations will no longer possess 
suitable energy invariants. 

However, in other studies, the thermal 
equation has been further simplified to be- 
come 


Tellus XII (1960), 4 


369 


where ©; is a standard value of ©, dependent 
upon p alone. Equation (23) contains no terms 
of class (1,3). If the expression for A is also 
simplified to become 


I 00 al + 
A-Laoast fr (229) Ve: Gran 


(24) 


a form resembling (4), it will follow from (23) 
that 


aA 
ns Hyp oe I p*-10,0,dM. 


(25) 

The rate of change of A is then identical 
with expression (22), the rate of change of 
P+I as determined from the unsimplified 
equation (15), so that the sum of kinetic 
energy and available potential energy will be 
conserved. For the purposes of many studies, 
this sum forms a sufficient energy invariant. 

However, equation (23) allows no variations 
of P+I, and hence no variations of S. If we 
wish to describe the static stabilization accom- 
panying the release of kinetic energy, and any 
consequent tendency to suppress the further 
growth of disturbances, we should retain 
all the terms of class (1,3) in the thermal 
equation (TS). 


4. Energy-preserving n-layer models 


In this section we shall establish a set of nu- 
merical prediction equations, for a model 
atmosphere in which the vertical dimension 
is replaced by a finite number of layers. We 
shall do this in such a way as to retain the 
various energy invariants. Accordingly, we 
may begin with one of the systems described 
in the last section. We shall use the equation 
of balance as a filtering approximation, and 
include the terms of class (2,3) in the vorticity 
equation. The further simplifications to be 
made if we wish to use the geostrophic equa- 
tion will be obvious. 

It will be convenient to 
variable 


introduce the 


(26) 


so that y: =2X/0p and w,=V 2X. If we omit 
the numerical .ubscripts, which are now super- 
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fluous, the thermal equation and the vorticity 
equation may be written. 


210) dX 0 
= Jy. 0)-9- (075°) 


oot = - J(y, f+ v2y) - v. (1 3 


OX OX 2) | 
-| v2wv + y? vy — v2Xv 
| D die re VU 
(28) 


while the equation of balance may be written 


gvtz=v-foy+o-| 2 ve (vp: |. 
(29) 


The reason for the particular grouping of 
terms in (27), (28), and (29) will soon be 
apparent. 

With the aid of (12), (13), and (14), the 
equation of balance may be converted into a 
generalized thermal wind equation 


EY ney à Loy aa de. à. 
CpPoo 9 d(p*) f y) 2 (p*) 
| vv = = v(vy: oy) (30) 


Equations (27), (28), and (30), together with 
the appropriate boundary conditions, form a 
closed system of three dependent variables ©, 
y, and x. 

The corresponding system with the geo- 
strophic equation, and without the terms of 
class (2, 3) in the vorticity equation, may be 
obtained simply by omitting the terms con- 
taining square brackets from (28) and (30). 

Let us now replace the three-dimensional 
atmosphere by n layers, bounded by the n+1 
isobaric surfaces to Pre = =F Spo numes 
bered from the ground upward. Thus py 
still represents surface pressure, while Pan = 0: 
The isobaric surfaces need not be spaced at 
equal intervals. Let odd subscripts from 1 to 
2n—1 denote the n layers. The mass of the 
atmosphere is now given by 


JdM=g LE" (pis = pyar) fdH, (31) 
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where I denotes a sum over all odd values of j, 
and dH is an element of horizontal area. 

We must now replace the system of di- 
fferential equations by a modified system in 
which finite differences replace derivatives 
with respect to p.-Our problem is to do this 
in such a way that reversible adiabatic proc- 
esses still have numerically equal effects upon 
kinetic energy, total potential energy, available 
potential energy, and gross static stability. To 
this end, we define © and y within each layer. 
At this point we depart from many of the 
currently used models in which the wind 
field is defined at n levels and the temperature 
field at n - 1 levels (see CHARNEY and PHILLIPS, 
1953). We define X at the surfaces separating 
the layers, so that in particular Xp =X, =O. 

The total potential energy and the kinetic 
energy are now given by 


Prlecp gt E (ps Pisa) pil Ada 
(32) 


and 


I : 
K== gt 2 (Bi-ı Bir) SI vydH (33) 


In order that (32) have meaning, however, 
we must have some rule, such as linear inter- 
polation, for defining p within the layers. 

The finite-difference forms of (27) and (28) 
may be obtained by replacing each indicated 
vertical derivative by a difference across a 
layer; thus 


99; XX 
18 I Oye yaaa 
Ot i» 9) te Piva pres 
ON es O;,:v?X; (34) 
Pi-ı — Pi+1 
Fa 
5 = Jp, f+ v2) + 


+ (pit Br) vf (Re 
+ (Bis Pro IV» [pv (X-ı - Xji1) + 
+ V9( Xj Kr) Vp 

a j-1 VPj-1 — VPX4 1 VYj41)| (3 5) 

This explains our grouping of terms in (2 
7) 
and (28); the vertical derivatives have been 
arranged so that X is referred to only at 
the surfaces separating the layers. However, 
in order that (34) and (35) have meaning, 
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we must have some rule, such as linear inter- 
polation, for defining © and y at the surfaces 
separating the layers. 

Upon integrating by parts, and again ob- 
serving that the divergence of any vector, and 
the Jacobian of any two scalars, vanish when 
integrated throughout the atmosphere, we 


find that 
2(P PIE 


Ste HS len 
ot fo Be UN pean) 


> fe +1 V2O;414H, 


(36) 


— =g 12 [Xrıv SY (Wy - Wire) d+ 
Ron SE AS 


[O7 VYj= V'Yj+a VYj+2) — 


I 
we v(vp- VY Via v.)| (37) 
provided that we let 


Wir => (y; + Yi+2) for odd j. (38) 


Comparing (36) and (37), we see that total 
energy is conserved provided that 


pee Ont (pepe live 
ae a) = (Pi Spee") atv: 


| (ow VYj— Visa VYj+) — 
if 
Lo (ey vw vus vw) | (35) 


Since this relation is a logical finite difference 
approximation to the generalized thermal wind 
equation (30), we have a set of equations with 
an energy invariant. 

From equation (34) it follows that, 


One 
— Z'(f-1-pi+1)0;dH=0, (40) 


dt 


so that the equations conserve the average 
value of ©, and hence conserve the difference 
between gross static stability and kinetic energy. 

In general, the equations cannot conserve the 
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entire statistical distribution of ©. Neverthe- 
less, since it follows from (34) that 


Dr à 
2 (Pj-ı - pi) J 9j dH = 


=> Xj41(20;41 = 9; — Oj2)- 

(0; ~ 0.2) 4H, Ua) 
the average value of ©2 will be conserved if 
we let 


9;4,= (9; + Oj+2) for odd j. (42) 


Thus, although the sum of kinetic energy and 
available potential energy, as originally de- 
fined, is not conserved, the sum of K and 
a modified form of A is conserved. This modi- 
fied form of A is the excess of P+J above the 
minimum value of P+I which could accom- 
pany any mass distribution with the same 
average values of © and OX. 

We still have some freedom of choice, since 
the rule for determining p within the layers has 
not been specified. For definiteness, let 


I 
Pi => (Pi-ı + pj+1) for odd j. (43) 
The system of equations (34), (35), and (39), 
together with the auxiliary definitions (38), 
(42) and (43), is now complete. 

Of special interest is the case where n =2 and 
Pa =po/2, which may be the simplest possible 
numerical prediction model with variable 
static stability. It is convenient to use as 
dependent variables the mean potential tem- 
perature © and the static stability o, the stream 
functions y and r for the mean wind and the 
wind shear, and the velocity potential x of the 
lower layer, so that 9,=© +0, 0, =0 - 5, 
Pa=V+T, Yi=y-t, and X,=pox/2. The 
governing equations (34), (35), and (39) then 
become 

20 


Er — J(y, 9) -J(t, o) +v-0vx, (44) 


D = Im) -J(e, 6) + 99-94, (49) 
I otp = ~J(y, Vy +f) Ja ver) + 


+v.[v2rvy + viyvrl (46) 
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Pra ily, HOT vy A) + 


+ v.fvg+v-[wypv2] (47) 
b,?O=v.fvr+Vv- 
[v?pvr+v?rvp- v(vy- vo), (48) 


where, because of (43), 


1) (Tree 0 


The corresponding system with the geo- 
strophic equation, and without the terms of 
class (2, 3) in the vorticity equation, is ob- 
tained by omitting the terms containing square 
brackets from (46), (47), and (48). In problems 
where f may be treated as a constant, (48) 
then simplifies to bc,O=ft, so that the term 
- J(r, ©) in (60) drops out. 

The total potential energy is given by 


P+I=potpg-1f(aO -bo)dH, (so) 


a -- [el 1. (2)] =0.797. (51) 


The kinetic energy is simply 


where 


K=<pog tf (vp: 9p + vr vo)dH. (52) 


The gross static stability should be a quantity 
dependent on 5, and obtainable by adding a 
multiple of © to - (P+T), where again a bar 
denotes a horizontal average. It is therefore 


given by 


S = bpoc,g"? fodH, (53) 


the negative of the second term in (so). 
Finally, the mean-square potential tempera- 
ture, given by 


As ar A? + SANTE ee aT ñ 
Gig=O@ +0 #0 7% 0,0. (54) 
is conserved, where 0’ = 9 - @ and 6 =o -5. 


Since © is also conserved, 6 has an absolute 
maximum 6», given by 


On = 0° +O? + 0%. 


(55) 
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The available potentiäl energy is then the 
excess of P+I, above the value of P+I ob- 
tained by substituting 9, for a in (so), i.e., 


A = bpotpg* J (Gm - 6) dH. (56) 
It then follows from (55) that 


_ [9% +02 
A= bpo%g Due. dH. (57) 


Hence, as in expression (3), A is given by a 
weighted average of the variance of potential 
temperature within isobaric surfaces. 

We thus have a simple two-layer model 
which properly describes the relations between 
total potential energy, kinetic energy, available 
potential energy, and gross static stability. 

Finally, we note that the model may be 
reduced to what is essentially one of the fami- 
liar two-layer models simply by discarding 
equation (45) for do/dt, replacing it by the 
relation o= constant. The latter model will 
preserve the sum of kinetic energy and 
available potential energy, but will not 
describe the relation between static stability 
and energy. 


6. Uses of the simplified equations 


During the past few years so many multi- 
layer models, and particularly two-layer mod- 
els, have been devised for numerical prediction 
that it might hardly seem worth while to add 
still another model to the collection. Indeed, 
the two-layer model presented in the previous 
section could probably not be justified on the 
grounds that it should yield better short-range 
forecasts, since the lack of variable static 
stability in other two-layer models is proba- 
bly not the primary reason for the errors in 
prediction. Such problems as improper side- 
boundary conditions and inadequate represen- 
tation of the initial three-dimensional wind 
and pressure fields are still present. 

The chief value of the model, then, is likely 
to be found in theoretical studies of the general 
circulation or similar circulations. For this 
purpose, additional terms should be appended 
to the equations, to represent the affects of 
heating and friction. | 

The two-layer model, with heating and 
friction, should be suitable for studying the 
flow in the “dishpan” experiments (cf. Furtz 
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1953). Here it would be relatively easy to 
solve the nonlinear equations for the steady 
symmetric flow in equilibrium with symmetric 
heat sources and sinks. Attempts to solve the 
Navier-Stokes equations for such a flow have 
been made by Davıss (1953) and others; great 
difficulties were encountered except when the 
equations were linearized. Once the symmetric 
flow is determined, it can presumably be 
tested for stability by the usual perturbation 
procedure. 

The two-layer model should also be suitable 


for studying various features of the general 


3713 


circulation, particularly those features which 
are also found in the dishpan. A special form 
of this model has already been used by Bryan 
(1959) to investigate some characteristics of 
the energy cycle. 

Problems involving the connection between 
the troposphere and the stratosphere might be 
studied with a three-layer model. The extent 
of the linkage between the troposphere and 
very high levels might be investigated with a 
model of several layers, which successively 
decrease in mass from a thick lowest layer to a 


thin highest layer. 
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The Exchange of Kinetic Energy between Larger Scales 
of Atmospheric Motion 


By BARRY SALTZMAN and AARON FLEISHER,! Massachusetts Institute of Technology 


(Manuscript received March 21, 1960) 


Abstract 


The rate of transfer of kinetic energy between the different scales of eddies of the soo mb 
geostrophic wind has been measured for each day of 1951. In the mean, there is a net transfer 
from intermediate (cyclone) scales to both longer and shorter scales. 


The rate, per unit area and per unit pressure from eddies of wave number n to eddies of 
difference, at which kinetic energy is transferred all other wave numbers is (SALTZMAN 1957) 


L*(n)=L*(n) + L"(n) (1) 


Po 7/2 © 


at ER ug {U (m) LUC à Un m) - UM) U(-n = ni] 


“th ewe 
+ V(m) [V(-#) Un - m) - V(n)U(-n- m)}} 
Pre gt n)[U (m) 7 (n-m) cos #]s + U(n)|U(m) V(-n-m) cos 6]; 
+ V(-n)[V(m)V(n-m) cos $]s + V(n) [V(m)V(-n-m) cos ]s} 
= 22-8 (y(n) [U( =) U(n -m) + U(n)U( = n =m) 


- U(m)[V(- n) U(n - m) + V(n) U(-n- Dh cos ddddp 


Po nf © 


{ut [Um One rem Tune 


+ V(-n[V (m) 2(n-m)],+V (n)[V(m)Q(-n- m)! 
cos Pdddp 
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A = longitude 
& = latitude 
p = pressure 


t =time 


ce) 


m, n = wave number around a latitude circle 


1 (A, d, p, t) = (a cos d) dajdt 
v (A, ,p, N 
o (A, #, p, t) = dp| 

F (x 


ab 


2% 


F(u) 
F(v) 
F 


Je 
= 
iam 
1) =F) 


ee 
sat 


= zonal wind component 
= meridional wind component 


= individual pressure change 


(2 ae x (A, 6, p, t)e-'"?dA = Fourier transform of x 


Po = pressure at the earth’s surface 


a=radius of the earth 


g = acceleration of gravity, 


(x)E=ax/OE. 


Since ZL*(n)=o we can think of L*(n) as 
zen 
the rate at which the existing kinetic energy 


is being redistributed among all the scales 
of atmospheric eddies, this process being in- 
dependent of other processes which can alter 
the total eddy kinetic energy—i. e., transfers to 
or from the zonal motion (n= o), viscous 
dissipation, and conversions to or from poten- 
tial and internal energy. 

We are concerned now with estimating 
L’(n), the contribution to L*(n) from inter- 
actions of the horizontal winds, for the scales 
represented on hemispheric charts. For this 
purpose we apply the geostrophic approxima- 
tion. We do not think these charts justify a 
resolution in more than 15 wave numbers and 
therefore take all scales of motion smaller 
than wave number 15 to be zero. We label 
the result of this approximation L(n). Then, 


X L(n)=o (2) 


which is to say that L(n) measures the geostro- 
phic redistribution of kinetic energy among 
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only the fifteen wave numbers. Transfers due 
to non-geostrophic components and transfers 
between the wave group # =1 to 15 and shorter 
waves are considered to be separate processes 
and are not treated here. 

Our data are the 500 mb northern hemi- 
sphere charts for every day of 1951. Contour 
heights are given at every ten degrees of 
longitude and every five degrees of latitude 
from 15° N to 80°N. These are the data from 
which we obtained the spectrum of the rate 
of transfer of kinetic energy between the 
zonal current and the eddies (SALTZMAN and 
FLEISHER, 1960). 

We assume that events at 500 mb are rep- 
resentative of the vertical average, at least 
in regard to sign. 

Our finite difference scheme limits the inte- 
gration to the region between 22.5° N and 
72.5° N. For these new boundaries, 


u al 


(brackets denote a zonal average and primes 
the deviation therefrom). However, if our 


b= 22.5 


b=72.5 
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results are to be representative of the globe it is 
necessary that (2) be satisfied, which means 
that the net eddy transport of eddy kinetic 
energy into this region must be zero. We 
have measured this and found it quite small. 
In place of a sufficient condition we assume 
that the region we treat is large enough to 
represent the interactions over the globe. 

The integral involves a summation of triple 
products of Fourier components—effectively a 
triple correlation. The variability of the 
atmosphere and the quality of the data hardly 
encourages such an enterprise. It should be 
expected, therefore, that only a long-time 
mean will have a variance sufficiently small to 
be worth quoting. We present the annual mean 
and, in computing its error, take half of the 
365 days to be independent. To obtain addi- 
tional smoothing, we consider three groups of 
wave numbers: ‘long’ waves I =n = 5, ‘cyc- 
lone’ waves 6 =n = 10, ‘short’ waves II Zn = 
=15, and we denote the transfer function 
for these groups by L;, i = 1, 2, 3, respectively; 
1.e., 


Ly= 2 L{(n) 
n=1 
10 
L,= 2 L(n) 
n=6 
15 
Les SL (n) 


jo at 


The loss in resolution is appreciable, but phys- 
ically less than might be supposed since the 
same wave number at different latitudes cor- 
responds to different physical scales. 

Results 


We define 
x = annual average of x 
o (x) = (x? — x*)} = standard deviation of x 


29 (x) = error of x 
V182 


The mean values obtained and their errors, 
INvTO À ergssecr cm /mbhedare 
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Fig. 1. Estimate of the mean rate of kinetic energy ex- 

change between long waves (n=I to 5), cyclone waves 

(n=6 to 10) and shorter waves (n=II to 15), based on 

daily s00 mb geostrophic measurements for the year 
1951. Units are 10-3 ergs sec + cm? mb. 
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Fig. 2. Estimate of the mean rate of kinetic energy ex- 
change between the zonal current -(n=o), long waves 
(n=1to 5), cyclone waves (n=6 to 10) and shorter waves 
(n = II to Is), based on daily soo mb geostrophic meas- 
urements for the year 1951. Units are 107 ergs sec-1 
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and 
AL; = + 50+ 56 


from which we conclude that L, and L, are 
significantly positive, L is significantly nega- 
tive, and as required, YL; is not significantly 
t=1, 2,3 
different from zero. In physical terms, there is, 
in the mean, a net transfer of kinetic energy 
from the ‘cyclone’ band (n=6 to 10) to the 
long-wave band (n=1 to 5) and to the short- 
wave band (n=ı11 to 15). Fjortort (1953) 
has demonstrated the possibility of transfers 
of this kind, and SALTZMAN (1959) has con- 
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structed a model which displays these energy 


flow characteristics. 
To estimate the numerical magnitude of 
this effect we assume that ZL; is entirely in 


t 
the nature of a computational noise which 
affects each of the three wave-groups equally. 


Thus, we subtract 1/3 ZL; froma Land J. 
to achieve the balance required by (2). The 


results are pictured in Figure 1. An arrow in- 
dicates a net loss or gain to or from the other 
two wave groups. 


In Figure 2 we display the results shown in 
Figure 1 combined with those presented by 
SALTZMAN and FLEISHER (1960) for the inter- 
actions between the mean zonal current and 


I 
the eddies. This is our estimate of the mean 
geostrophic energy flow in the domain n 
=0to 15. Estimates of the release of potential 
and internal energy in this domain are given 
by SALTZMAN and FLEISHER (1959) and Wan- 
NIELSEN (1959). There are no direct measure- 
ments of the transfer of kinetic energy between 
these scales of disturbances and scales of wave 
number larger than rs. 

We are indebted to the people of the Com- 
putation Center of the Massachusetts Institute 
of Technology for time on the IBM 704 and 
for running our program. This research has 
been sponsored by the Geophysics Research 
Directorate of the Air Force Cambridge Re- 
search Center under Contract No. AF 19(604)- 
2242. 
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Abstract 


Multiple regression equations for predicting s-day mean surface temperature in the United 
States from prior anomalies of 700-mb height and surface temperature are derived by a screening 
procedure for fall, winter, and spring. The heights required for an objective forecast are approxi- 
mated with the aid of daily numerical prognoses prepared at the soo-mb level by a barotropic 
model and at the 700-mb level by a baroclinic model; the required temperatures are estimated 
by combining observed daily values with short-range prognoses. The resulting objective tem- 
perature, predictions are verified on independent data at a number of cities, and the sources of 
error are analyzed. Comparison is made with forecasts prepared by local temperature per- 
sistence, continuity of error, and conventional methods. It is concluded that the objective 
method can produce in a short time temperature predictions which have appreciable skill 
beyond persistence. Thus a valuable aid and base for further improvement has been made 


available to the forecaster. 


1. Introduction 


Considerable research on the relation between 
temperature and circulation for extended pe- 
riods was conducted in the Extended Forecast 
Section of the United States Weather Bureau 
about ten years ago by Martin and LEIGHT 
(1949) and Martin and Hawkins (1950). 
During the past few years interest in this 
problem has been revived because of the 
availability of high-speed electronic computers 
(Kettey and MALONE, 1955). In a previous 
paper (KLeIN, Lewis, and ENGER, 1959) a sta- 
tistical screening procedure was used to obtain 
multiple regression equations ‚for predicting 
5-day mean surface temperature as a linear 
function of selected 700-mb heights. Applica- 
tion of these equations to heights obtained from 
barotropic prognoses produced temperature 
predictions of positive skill during the test 


winter of 1957—58. The forecasts were con- 
siderably improved by including as a predictor 
the local value of current 5-day mean surface 
temperature, 

The last paper described above was primarily 
exploratory in nature and deficient in several 
respects. The only temperature predictor used 
in the regression equations was the local tem- 
perature; i.e., at the station for which the fore- 
cast was made. Furthermore, the final equations 
were derived on desk calculators without use 
of the screening technique and were therefore 
limited to 11 cities, 1 season, and 3 predictors. 
The resulting objective forecasts were tested 
on only 16 independent cases. 

The present paper will describe an objective 
method which did not suffer from the afore- 
said limitations. In addition to persistence at 
the local station, temperatures at a large 


1Part of this paper was presented on May 7, 1959 at the joint meeting of the American Meteorological 
Society and the American Geophysical Union in Washington, D.C. i 
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network of cities in the United States were 
considered as possible predictors. By use of a 
high speed electronic computer, the screening 
technique was applied to the field of surface 
temperatures as well as 700-mb heights. Re- 
gression equations were derived for 3 seasons 
and 39 cities, and as many as 7 predictors were 
included in each equation. The forecasts were 
tested on a large sample of independent data 
consisting of a total of 226 cases. 

We shall first discuss the objective method, 
including the data used to derive it (section 2), 
its meteorological interpretation (section 3), and 
the screening procedure (section 4). A series 
of experiments with various numerical prog- 
noses as input data will be described in section 
5.The results of applying the objective method 
during a 2-year trial period will be presented 
in sections 6 and 7, and some comparison will 
be made with forecasts produced by other 
methods. Limitations of the present objective 
method and suggestions for improvement 
will be discussed in section 8. The principal 
findings are summarized in section 9. 


Fig. 1. Network of 70 grid points used to delineate the 

700-mb height field. Heights were taken at standard 

intersections of latitude and longitude, 10 degrees apart, 
marked by open circles. 


2. Data Used 


There were two basic parameters used in 
the forecast scheme. The first was the departure 
from normal of 5-day mean 700-mb height 
over North America and adjacent oceanic 
areas. Heights? were taken at a network of 70 
grid points, 10 degrees apart, extending from 
30° to 70° north and from 50° to 180° west 


(figure 1). 


2 All references to heights and temperatures throughout 
this paper should be understood to apply to their anom- 
alies, or deviations from local normal, rather than to 
their absolute values. 
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Fig. 2. Network of 39 stations used for surface tempera- 
ture. 


The second parameter was the anomaly of 
5-day mean surface temperature at a network 
of cities covering all parts of the country. For 
fall and spring the 39 stations shown in figure 2 
were used, but for winter only 30 were avail- 
able (see fig. 2 of KLEIN and others, 1959). In 
selecting cities, consideration was given to 
such factors as homogeneity of record, repre- 
sentativeness of exposure, and availability of 
teletype data. 

In order to clarify the time periods involved 
in the forecasts, a schematic calendar has been 
prepared, as shown in figure 3. Let us call 
the day on which the prediction is made day o. 
The 5-day mean forecasts routinely prepared 
in the Extended Forecast Section of the United 
States Weather Bureau apply to the period 
beginning 2 days after forecast day (day 2) and 
ending 6 days later (day 6); i.e., to the s-days 
centered 4 days after forecast day and there- 
fore designated T',. Since it is difficult to pre- 
dict accurately the 700-mb height field for a 
period this far in advance, the 5-day mean 
heights required for the objective forecast were 
taken from a period 2 days earlier, extending 
from day o to day 4. These heights are desig- 
nated H, because they are centered 2 days 
after forecast day. This period was chosen 
because it gave good results in our earlier 
study (KLEIN and others, 1959), when it was 
noted that surface temperatures tend to lag 
behind upper level heights, due in part to 
the time required for advection to operate. 

In operational practice, the daily barotropic 
forecasts prepared by the Joint Numerical 
Weather Prediction Unit (JNWP) at Suit- 
land, Maryland (CRESSMAN, 1959) were com- 
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DAY 
—2 


= 


FORECAST DAY —— O 


5 


6 


Fig. 3. Schematic calendar with various s-day mean 
periods depicted relative to forecast day. 


bined into a s-day mean map by simply 
averaging the last observed map (day 0) with 
the forecasts for each of the next 4 days from 
day 1 to day 4. The resulting objective estimate 
of H, was readily available because it has been 
prepared on a routine basis in the Extended 
Forecast Section for the past few years, as ex- 
plained by Namras and Collaborators (1958). 


The temperature predictor used in the 
objective method was centered on forecast 
day and therefore designated Ty. It covered 
the period from 2 days prior to forecast day 
to 2 days later and thus overlapped T, by one 
day (day 2). Half the data required for T', had 
already been observed, and the other half was 
estimated from forecasts for the next 2% 
days routinely prepared at Weather Bureau 
district forecast centers throughout the United 


lized s-day mean surface temperatures centered 
on forecast day and 700-mb heights centered 
2 days later to predict 5-day mean surface 
temperatures centered 4 days later, or: 


Th =f (To, H;) (1) 


In the forecast routine, H, and T, were 
estimated from the numerical and district 
forecasts. Since these estimates were not 
available for a long period of record, the 
research was conducted with only observed 
values of both quantities. The fall season was 
defined as the 14-week period from mid- 
September to mid-December, and one case 
was chosen each week over a 10-year period 
from 1948 to 1957, making a total of 140 cases. 
A similar procedure was followed for the 
other seasons, so that 140 winter cases were 
selected from mid-December 1947 to mid- 
March 1957, and 140 spring cases from mid- 
March 1949 to mid-June 1958.3 


3. Synoptic Implications 


The first step in deriving the objective meth- 
od was to obtain simple linear correlation 
coefficients between 5-day mean surface tem- 
perature (Ty) at each city plotted in Figure 2 
and s-day mean 700-mb height 2 days earlier 
(H,) at each of the 70 grid points shown in 
Fig. 1. The resulting correlation fields are illus- 
trated in figure 4 for New York City during 


3’ A more refined treatment would have consisted of 
standardizing the data, so as to equalize the variability 
of different months within the same season. Better yet 
would have been to subdivide the data by months, 
rather than seasons, but an insufficient number of years 
was available for this purpose. 


Wo 
NR 


SHE 


Fig. 4. Simple linear correlation coefficients between 

s-day mean surface temperature anomaly at New York 

City (located by asterisk) and s-day mean 700-mb 
height anomaly 2 days earlier for 140 fall cases. 
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the fall season. The isopleths of equal correlation 
indicate that fall temperatures at New York 
during the 10 years studied were most closely 
related to 700-mb heights in the immediate 
vicinity. The magnitude of the correlation was 
-60 and the sign was positive, thereby showing 
that above normal heights were associated 
with above normal temperatures, while low 
heights were followed by cold weather. This is 
a reflection of the fact that the troposphere 
usually approaches an equivalent-barotropic 
state (CHARNEY, 1949), with the contours tend- 
ing to be in phase with the thickness lines. 
In addition, surface temperature is generally 
well correlated with thickness in the lower 
troposphere (ROWE, 1944). 

There is another important center of maxi- 
mum correlation in figure 4. In the Northern 
Rocky Mountain States the highest correlation 
was .53 but the sign was negative, thus indi- 
cating that above normal heights in this area 
were followed by below normal temperatures 
at New York, while low heights accompanied 
warm weather. One synoptic interpretation 
is as follows: Since the Rockies are normally 
about half a wave length upstream from New 
York, a ridge of large amplitude in the West 
usually results in strong northwesterly flow 
of cold polar continental air into New York; 
conversely, low heights over the Rockies imply 
westerly or southerly flow in the East and thus 
milder air masses. 

Figure 5 Asummarizes the results obtained by 
constructing similar correlation fields for all 
39 cities during the fall season. The heavy dots 
show the points of maximum correlation 
between 700-mb heights and surface tempera- 
tures 2 days later at cities located near the tip 
of the arrows. For example, at Bismarck, 
North Dakota, the maximum correlation was 
—.64 with 700-mb heights in the Alaskan 
peninsula. Similarly, most cities in central 
and northwestern sections of the United States 
had points of highest correlation about 1,000 
to 2,000 miles to the northwest. On the other 
hand, most cities in the East and Southwest 
were like New York, with the maximum 
correlation positive and located close to the 
station. All cities were basically similar, 
however, in exhibiting one center of positive 
correlation near the station and one center of 
negative correlation about half a planetary 
wave length upstream. 
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Fig. 5. Maximum simple linear correlation between 5-day 

mean anomalies of surface temperature and 700-mb 

height 2 days earlier for A) fall, B) winter, C) spring. 

The point of highest correlation is shown as a heavy 

solid circle, the city is located at the tip of each arrow, 

and the value of the correlation is written next to each 
city. 


The above statements apply equally well to 
the winter season, for which the corresponding 
diagram is presented in figure 5 B. The marked 
similarity between figures 5A and 5 B suggests 
that fall and winter tend to be climatologically 
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Fig. 6. Simple linear correlation coefficients between 

s-day mean surface temperature anomaly at New York 

City (starred) and s-day mean surface temperature 

anomaly at network of 39 cities 4 days earlier for 140 
fall cases. 


homogeneous with respect to the interrelation 
between temperature and circulation. In spring, 
on the other hand, a somewhat different rela- 
tionship prevails since at all but three cities the 
correlations in figure $C are positive, and the 
arrows originate in the vicinity of the station. 
This indicates an increased importance of local 
height, at the expense of heights upstream, 
during the spring months. This characteristic 
is even more pronounced in summer, according 
to an earlier study by Martin and HAWKINS 
(1950), when local factors like cloudiness and 
vertical stability have a greater effect on surface 
temperature than does advection of air masses. 

A second step in deriving the objective 
method was to obtain simple correlations 
between s-day mean temperature at each city 
(T,) and s-day mean temperature for a period 
(overlapping by one day) centered 4 days 
earlier (To) at other cities. For example, the 
correlation field of figure 6 shows that at 
New York City during the fall the persistence 
of local s-day mean temperatures 4 days apart 
gave a correlation of only +.33. Higher 
correlations were obtained by considering the 
temperature 4 days earlier at points to the 
west, and the maximum correlation (+ .54) 
was given by a point in the upper Great Lakes. 
The positive sign shows that low temperatures 
in the Great Lakes tend to be followed by 
cold weather in New York 4 days later, and 
conversely for warm weather. This is an 
expression of the normal westerly steering 
flow prevailing at New York and may be 
associated with the fact that s-day mean 


troughs and ridges move slowly eastward on 
the average (NAMIAS, 1947). 

Figure 7A summarizes the results of a simi- 
lar analysis at all 39 cities during the fall 
season by showing the points of maximum 
correlation between temperatures 4 days 
apart. For example, at Albany, New York, 
the results are similar to those at New York 
City, with the highest correlation being + .57 
with temperatures in the Upper Lakes region. 
The distance to the points of maximum corre- 
lation ranged from a few miles for cities in 
the Great Basin to as much as 800 miles for 
eastern cities like Jacksonville, Florida, and 
Buffalo, New York. It is interesting to note 
that many of the arrows point in a southerly 
as well as an easterly direction, suggesting 
that 5-day mean temperature anomalies tend 
to move southeastward over intervals of 4 
days. This result may be due in large part to 
the strong tendency for cold polar anticyclones 
to move southward as well as eastward in 
North America (KLEIN, 1957). Figures 7 B and 
7C show that similar relations apply during the 
winter and spring seasons. 


4. Screening Procedure 


Thus far we have considered simple linear 
correlation coefficients which tell how well 
we could forecast temperature by using only 
one predictor. The next task is to improve 
these results by selecting additional predictors 
which will contribute significantly and inde- 
pendently to the temperature. This was done 
by a step-wise method of multiple regression 
known as the screening procedure which had 
been programmed for the IBM-704 electronic 
computer by Mixer (1958). As explained in 
our previous paper (KLEIN and others, 1959), 
the machine automatically selected the second 
most important variable by maximizing the 
multiple correlation between temperature 
and the combination of the first predictor plus 
any other. This process was repeated, and 
additional variables were selected one by one. 


For example, at New York City in fall the 
most important single variable was the 700-mb 
height 2 days earlier at the intersection of 40° N 
and 70° W. The predictor which contributed 
the most additional information was the 
surface temperature 4 days earlier at Interna- 
tional Falls, Minnesota. Combination of this 
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Fig. 7. A) Maximum simple 
linear correlation between 
s-day mean anomalies of 
surface temperature 4 days 
apart for 140 fall cases. The 
point of highest correlation 
is shown as a heavy solid 
circle, the city is located at 
the tip of each arrow, and 
the value of the correlation 
is written next to each city. B) 
Same as A) for 140 winter 
cases. C) Similar to A) for 
140 spring cases, except that 
the data were obtained from 
the first station temperature 
selected by the screening pro- 
gram rather than from the 
analyzed correlation fields. 
The long arrows originate 
at the predictor city and ter- 
minate at the predictand city, 
the short arrows indicate that 
the city itself was selected, 
and the squares indicate that 
no station temperature was 
selected as a predictor. 
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variable with the first yielded a multiple 
correlation of .68. The maximum multiple 
correlation using the combination of these 
first two predictors and any other one was 
.75, produced by adding the 700-mb height 2 
days earlier at 50° N, 110° W. Since no addi- 
tional variable was able to increase the ex- 
plained variance by more than 2% %, the 
screening process was stopped with these 
three predictors, all of which were located 
near the centers of maximum on the original 
correlation fields (figure 4 and 6). The linear 
multiple regression equation used in preparing 
the objective forecast during the fall season 
can then be written as: 


T, (New York) = .75 + 
+.104 H,(40 N, 70 W) +.263 Ts 
(Int. Falls) —.097 H, (50 N, 110 W) (2) 


where temperatures are in degrees Fahrenheit 
and heights in tens of feet. 


Multiple regression equations of the type 
illustrated in equation (2) were derived sepa- 
rately for fall, winter, and spring at each of the 
cities shown in figure 7. Physical interpretation 
of these equations was similar to the discussion 
already given for New York. In nearly all 
cases at least 2 heights were selected, one 
located in the vicinity of the remote center of 
correlation to the northwest, and the other 
near the local center. In addition, at least 
one temperature was usually selected, generally 
from the center of maximum correlation near 
or northwest of the city. At most cities the 
screening procedure was continued after the 
third step, with one to three additional heights 
and temperatures (generally to the west of 
the station) being selected as significant predic- 
tors on the basis of synoptic reasoning as well 
as statistical tests. The final equations con- 
tained from 2 to 7 variables (heights and tem- 
peratures in about equal proportion) and ex- 
plained an average (for all cities and seasons) 


of 60 percent of the temperature variance of 
the developmental sample. 


5. Tests on Independent Data 


Since standard statistical tests of significance 
are not strictly applicable to this type of work, 
it is important to test the results on independent 
data. This was done by programming the 
multiple regiession equations for the IBM- 
704 and then making temperature forecasts 
at each of the 39 cities on each Thursday for 
12 weeks of the 1958 fall, from October 2 to 
December 18, for a total of 468 forecasts. The 
forecasts at each city were verified in terms of 
the explained variance. 


x (F; = 0;)? 
EV=1- — (3) 
E(N-0)% 


where F; is the forecast, o; the observed, and 
N the normal temperature. Values of EV 
were averaged for all 39 cities; the results 
are shown in table 1. 

The first test was made by using only 
observed s-day mean values of 700-mb height 
2 days earlier (H,) and surface temperature 4 
days earlier (T,) as predictors of T,. Since 
this was exactly the same type of data used 
in the original derivation, it provided a test 
of the stability of the regression equations. The 
average EV obtained was 54.7 percent (line 1, 
table 1), which is almost as good as the EV of 
60 percent explained on the original sample 
of 10 years (section 4). This result indicates 
that the regression equations are stable features 
which hold up well on new data. 

Thus far we have discussed observed values 
available after the fact. How would the equa- 
tions fare under actual operating conditions; 
i.e., using only estimated values of 5-day mean 
height and temperature available on forecast 
day? The required values of H, and Ty were 


Table 1. Percent of Variance of 5-day Mean Temperature Explained by Multiple Regression Equa- 
tions during Fall of 1958 (mean of 12 cases at 39 Cities). 


ap © D A 


Observed 700-mb height, observed surface temperature 
Barotropic 500-mb height, estimated surface temperature 
Observed 700-mb height, estimated surface temperature 
Observed 500-mb height, estimated surface temperature 
Baroclinic 700-mb height, estimated surface temperature 
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obtained with the aid of the barotropic height 
(0000 GMT) and district temperature forecasts 
as explained in section 2. The second line of 
table 1 shows that only 9.3 percent of the tem- 
perature variance at the 39 cities was explained 
in this way, a decrease of 45 percent from the 
results obtained in the first test with observed 
heights and temperatures. 

In order to investigate the sources of this 
decrease, two additional tests were made. In 
the first, estimated values of temperature were 
substituted for observed values, but observed 
7oo-mb heights were retained. The percent of 
variance explained was 49.2 (line 3), only 5 
percent less than the EV obtained in the first 
test using observed temperatures. This result 
indicated that the estimated temperatures were 
good approximations to the observed tempera- 
tures for the five days centered on forecast 
day (To). It therefore appears that most of the 
errors in the temperature forecasts were due 
to errors in estimating the heights. A similar 
conclusion can be drawn from the results 
obtained with a different sample in our earlier 
study (KLEIN and others, 1959). 

One source of height error is the fact that 
the barotropic forecasts were made at the 
soo-mb level, while the objective method was 
derived from heights at 700 mb. In all cases 
the required conversion was made by simply 
multiplying the soo-mb height by the ratio 
of normal 700 to soo-mb height for each 
month at each point.* To evaluate the inaccu- 
racies of this conversion, a second test was 
made by using 5-day mean heights observed 
at the soo-mb level and reduced to the 700-mb 
level, in conjunction with the estimated tem- 
peratures. Only 36.4 percent of the variance 
was explained in this way (line 4), compared 
to 49.2 percent for observed 700-mb heights. 
The difference of 13 % can be attributed to 
improper reduction from 500 to 700 mb. 


In an attempt to diminish this error, the 
regression equations were applied to baroclinic 
forecasts of 700-mb height rather than to 
barotropic forecasts of soo-mb height. During 
1958 JNWP introduced into its operations a 
two-level model which predicts the thickness 
between 850 and soo mb by taking account 


4230 symbols: 
N; 
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of horizontal advection and vertical motion. 
This prognostic thickness is subtracted from 
the barotropic soo-mb height to yield a prog- 
nostic 85o-mb height, and the 700-mb height 
is then obtained by interpolation. The equation 
currently used by JNWP is: 


LAS 37 Lg 70825 (4) 


where Z is the height in feet at the 700-mb, 
soo-mb, and 850-mb levels, respectively. 


Since the baroclinic model is run routinely 
up to only 36 hours in advance, it was not 
possible to obtain directly a 5-day mean map 
for the H, period at the 700-mb level. However, 
a good approximation could be obtained from 
the 700-mb 36-hour prognostic chart routinely 
prepared by JNWP at 1200 GMT on forecast 
day since it is valid on day 2 (fig. 3). This is 
the middle day of the required 5-day period, 
and therefore its map resembles the 5-day 
mean more closely than any other daily map 
(NAMIAS, 1947). Because the variability of daily 
maps is greater than that of 5-day mean maps, 
a damping factor of .5 was applied to the prog- 
nostic daily height anomalies before they were 
put into the regression equations. This is equiv- 
alent to taking H, as the mean of the 36-hour 
prognosis and the normal height at each point. 
Use of this H, as input data, again in conjunc- 
tion with the estimated mean temperatures, 
led to a marked improvement. An average of 
21.0 percent of the temperature variance was 
now explained (line 5), more than double the 
9.3 percent explained by the barotropic 500-mb 
maps. 


Part of this improvement was due to the 
fact that it was operationally feasible to base 
the 700-mb prognostic chart on observations 
made 12 hours later than those used as the 
starting point of the soo-mb mean map. In 
order to analyze additional causes of the im- 
provement, experiments were performed with 
daily barotropic soo-mb prognoses prepared 
from 1200 GMT data on forecast day. Each 
of these was used as input for the regression 
equations (in conjunction with the same 
estimated value of T, used previously) to see 
how much of the temperature variance they 
could explain. The results are shown in figure 8 
for the initial observed map (time o) and the 
prognostic charts valid 24, 36, 48, and 72 
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PERCENT OF TEMPERATURE VARIANCE EXPLAINED BY NUMERICAL PROGNOSES 
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Fig. 8. Percent of total variance of s-day mean surface 
temperature anomaly explained by barotropic forecasts 
of soo-mb height from o to 4 days in advance. The curve 
is extrapolated (dashed) after 72 hours. The open circle 
represents the variance explained by the 36-hour numeri- 
cal forecast of 700-mb height. All values are averages of 
39 cities during 12 fall cases from October 2 to December 
18, 1958. 


hours later. The points are easily fitted by a 
smooth curve which has a maximum explained 
variance of 14.2 percent at 48 hours, followed 
by a rapid decline to (extrapolated) zero around 
96 hours. Thus, none of the daily barotropic 
forecasts was able to explain as much of the 
surface temperature variance as the 21.0 per- 
cent explained by the baroclinic 36-hour 700- 
mb map, which is plotted as an open circle in 
figure 8. 


6. Verification in Terms of Temperature 
Classes 


Guidance forecasts of 5-day mean tempera- 
ture issued to Weather Bureau District centers 
by official forecasters of the Extended Forecast 
Section are expressed in five classes; much 
below, near, above, and much above normal. 
The center three have a climatological proba- 
bility of 1/, each; the first and last, 1/, each. 
In order to express the objective temperature 
forecasts in terms of these classes, they were 
first multiplied by the reciprocal of the mul- 
tiple correlation at each city since: 


0, 
R= 
09 


(s) 


OBJECTIVE FORECAST 
NOV 1418. 1959 


A 


We 
3 Ra, ! 


RT 


OFFICIAL FORECAST 
NOV. 14 18, 1959 


MEAN TEMPERATURE 
DEPARTURE OBSERVED. 
NOV. 14.18, 1959 c 


Fig. 9. Forecast and observed temperature anomalies for 
the s-day period from November 14 to November 18, 
1959. Analysis is in terms of the five temperature classes 
customarily used in the U.S. Weather Bureau (A) is 
the objective prediction, (B) is the official forecast, and 
(C) is the verifying map. The objective forecast had a 
skill score of 28 with 97 % correct within one class; 
the official forecast had a skill score of 53 with 96 % 
correct within one class. 


where R is the multiple correlation coefficient, 
oy is the standard deviation of the forecast 
temperatures, and oy is the standard deviation 
of the observed temperatures (EZEKIEL, 1941). 
The regression forecasts were inflated in this 
way in order to equate their variability with 
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Table 2. Verification of 5-day mean forecasts in terms of 5 temperature classes at 100 cities in 
United States during 3 seasons (1958—1959).! 
Ne at Ne ee eal Ee NTS Se Te ee 


Forecast Method 


Skill Score 
Fall Winter Spring All 


Percent within 1 
class 
Fall Winter Spring All 


Percent in correct class 
Fall Winter Spring All 


———— © 


11,0,33:10.32.0W28:0 37.48 |179:92.70.0171:207443 
14.7 |33-.3 34.5 34.0 33.9 |78.8 79.3 77.8 78.6 
012000 48277628 30075 1078 0007187114 
11.9 | 34.7 33.6 27.3 31.7 | 77.2 75.9 69.2 73.9 
14/7285 705100353053 01008870 0877117102 
14.5 |34.3 33.1 32.0 33.1 | 765 76.2 74.7 75.8 
0.0 21.9 59-4 


a Objective-barotropic 

500-mb mean Tara 2100710 
2, Objective-baroclinic 

7oo-mb daily 14.4 14.9 14.6 
33 Persistence-estimated 

mean temperature 9:97 79137 76:8 
4. Continuity of 

error-graphical TO HET 7103 
5. Official-Extended 

Forecast Section 15.286153 1909 
6. Official — 1952—1957 T5 ON DETNIS O 
7 Chance-random guess 
1 Fall —32 Cases from September 30, 1958 to December 11, 1958 


Winter— 39 Cases from December 14, 1958 to March 12, 1959 
Spring — 39 Cases from March 15, 1959 to June 11, 1959 


that observed and thus prevent an excessive 
number of forecasts in the near normal cate- 
gory. 

Objective forecasts were prepared routinely 
three times a week on the IBM-704 from 
estimated values of T, and Hy; i.e., only in- 
formation actually available on forecast day 
was used. Two different estimates of H, were 
tested: the s-day mean maps prepared from 
the barotropic soo-mb forecasts at oo00 GMT, 
and the daily 36-hour forecasts prepared at 
the 700-mb level from the baroclinic model 
at 1200 GMT. The inflated forecasts (equation 
5) were converted into one of the five tempera- 
ture classes according to the appropriate class 
limits at each city (Namias, 1947). After the 
forecast class for each city was plotted on a 
map, lines were drawn by interpolation 
delineating the distribution of the five classes 
over the entire country. A sample of a good 
objective forecast is shown in fig. 9 A. 110 fore- 
cast maps covering three seasons were verified 
at a grid of 100 points fairly evenly spaced 
over the United States. 


The verification results are summarized in 
table 2 on a seasonal and over-all basis. The 
first line shows that when the regression equa- 
tions were applied to the barotropic soo-mb 
mean maps, an average (for all three seasons) 
of 31.4 % of the objective forecasts were in 
exactly the correct temperature class, for a 
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skill score® (BRIER and ALLEN, 1951) of 11.6. 
Credit is sometimes given for a near hit, such 
as forecasting much below normal when 
below verifies. On this basis, the last column 
shows that 74.3 % of the objective forecasts 
were either exactly correct or correct within 
one class for all seasons combined. Line 2 
shows that application of the objective method 
to the baroclinic 700-mb 36-hour forecasts 
resulted in an improvement to 33.9 % in the 
correct class, an average skill score of 14.7, 
and 78.6 % correct within one class. The 
scores in line 2 were not only higher than 
those in line ı but also more consistent, with 
little variation from one season to the next. 
The remainder of the table contains various 
controls against which to measure the accuracy 
of the objective forecasts. Line 3 shows the 
scores obtained by local persistence; i.e., using 
the estimated mean temperature (To) as a 
forecast of Ty. Although these scores were 
definitely higher than those expected by chance 
(line 7), they were consistently lower than the 
scores given by either of the objective fore- 


5 The skill score is 


R—E 
s-(7= 5} 100 


where R is the number of correct forecasts, T is the 
total number of forecasts, and E is the number expected 
correct by chance. S is 100 for perfect forecasts, and it is 
positive, zero, or negative for forecasts better, equal, or 
worse than chance forecasts. 
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casts. Thus the skill of the latter cannot be 
ascribed solely to persistence of a short-range 
forecast. 

Line 4 gives the scores obtained by a meth- 
od based upon simple two-parameter graphs 
(MARTIN and Hawkins, 1950) and corrected 
by applying continuity of error in the manner 
described by Namras and Collaborators (1958). 
This method was approximately equal in skill 
to the objective forecasts based upon the soo- 
mb barotropic mean (line 1), but was not as 
good as the objective method applied to the 
baroclinic 700-mb maps (line 2). 

Line 5 contains the verification of the official 
guidance forecasts issued by experienced fore- 
casters of the Extended Forecast Section. 
During the 9-month test period, these fore- 
casters were able to make use of objective 
temperature predictions based on the baro- 
tropic soo-mb mean maps, local persistence, 
and continuity of error; but the objectives 
made from the baroclinic 7oo-mb 36-hour 
charts were not available to them. As a result, 
the scores in line 5 were decidedly superior 
to those in lines 1, 3, or 4 but merely equal to 
those in line 2. It is encouraging to note that 
both the official and baroclinic objective fore- 
casts were superior to the average of the 
official forecasts during the 5 years from 1952 
to 1957 (line 6), especially when near hits 
were counted (last column). 


7. More Recent Results 


Numerous attempts have been made to 
improve the results described in the preceding 


two sections by obtaining better forecasts of 
the 700-mb mean circulation. Although the 
details of these experiments will not be pre- 
sented here, they can be summarized by stating 
that in no case has any significant increase in 
accuracy been attained. The results clearly indi- 
cate that the 36-hour forecast of 700-mb 
heights made by JNWP at 1200 GMT of 
forecast day is. the best numerical input cur- 
rently available for the objective temperature 
method described in this paper. 

Beginning on September 15, 1959, objective 
temperature forecasts have been prepared on an 
operational basis, using these 700-mb heights, 
instead of the soo-mb mean barotropic heights, 
as input to the regression equations. The result- 
ing predictions have been verified in terms o 
temperature classes (section 6) for 116 cases 
through June, 1960. The results are summarized 
in table 3, along with comparative verification 
for other methods. 


Comparisons of tables 2 and 3 (line 2) shows 
that the objective temperature forecasts during 
the fall of 1959 were definitely superior to 
corresponding forecasts made the previous 
fall. This improvement was accomplished 
despite a sharp drop in the skill of the forecasts 
based upon either persistence (line 3) or con- 
tinuity of error (line 4). The objective predic- 
tions continued to exhibit superiority over 
both persistence and continuity of error during 
the winter and spring (table 3), despite a 
marked difference in the character of these 
seasons. The winter of 1959—60 was unusually 
persistent in nature, and all temperature fore- 


Table 3. Verification of 5-day mean forecasts in terms of 5 temperature classes at 100 cities in 
United States during 3 seasons (1959—60).! 


Skill Score 


Percent in correct class} Percent within ı class 


Line Forecast Method : mite : = 
eae gee Fall Winter Spring All| Fall Winter Spring All| Fall Winter Spring All 

L. Objective-barotropic 

500-mb mean Not computed 
By Objective-baroclinic 

700-mb daily 29,0 27.3 13.9 16.11 37.0 38.1 32.8 36.01,.87.87.83.,5, 180.87 82:0 
Br Persistence-estimated 

mean temp. 8.9 19.4 9.7 12.8) 29.3 36.8 29.6 32.0] 72.4 74.7 70.5 725 
4. Continuity of error- 

graphical 10.3 (8.0 10.7 13,416 31.125 35.062308, oe: an. 70. ; 5 
5. Official-Extended A Ne 

Forecast Section 20.4 27.0 17.2 21.6] 38.1 42.6 35.4 38.8| 79.6 85.3 79.3 81.4 
‘ Fall — 37 Cases from September 15, 1959 to December 8, 1959 


Winter — 40 Cases from December 10, 1959 to March ro, 1960 
Spring — 39 Cases from March 13, 1960 to June 9, 1960 
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OBJECTIVE FORECAST 
NOV. 19-23, 1959 


OBJECTIVE FORECAST 
NOV. 26-30, 1959 


OFFICIAL FORECAST 
NOV 19-23 1959 


OFFICIAL FORECAST 
NOV 26-30, 1959 


MEAN [TEMPERATURE 
OEPARTURE OBSERVED 
NOV. 19-23, 1959 


MEAN TEMPERATURE 
DEPARTURE OBSERVED 


Fig. 10. Objective forecast (A), official forecast (B), 

and observed (C) temperature anomaly classes for the 

s-day period November 19—23, 1959. The skill scores 

were 20 for the objective and 33 for the official, 

but both forecasts had 83 percent correct within 
one class. 


casts experienced a sizeable increase in skill. 
The spring of 1960 was rather difficult and 
changeable, and the accuracy of the forecasts 
diminished sharply. For all three seasons com- 
bined, the objective temperature predictions 
had a skill score of 18.1, a decided improve- 
ment over the corresponding score of 14.7 
obtained during the preceding year (table 2). 


Because many factors could have contributed, 
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NOV. 26-30, 1959 


Fig. 11. Objective forecast (A), official forecast (B) ‚and 

observed (C) temperature anomaly classes for the s-day 

period, November 26—30, 1959. The skill scores were 4 

for the objective and 14 for the official, the percentages 

correct within one class were 67 for the objective and 
75 for the official. 


and because only two years are involved, it is 
difficult to know to what extent this improve- 
ment will continue to hold in the future. 
Line 5 contains the verification of the guid- 
ance forecasts issued to the districts by official 
forecasters of the Extended Forecast Section. 
The average skill score of 21.6 was definitely 
better than the average skill score of the objec- 
tive forecasts (18.1), and this superiority was 
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maintained during each season. It was also a 
marked increase over the official score of the 
previous year (14.7, table 2), and the highest 
score on record for any 9-month period since 
1941. 

It is difficult to ascertain the exact causes of 
this increased skill over previous years. The 
official forecasts are based primarily upon 
prognostic charts of the 5-day mean 700-mb 
circulation for the H, period, plus numerous 
additional tools, as explained by NaMras in 
several papers (1947, 1958). One of these tools 
is the objective temperature forecast based on 
the 36-hour 700-mb baroclinic heights. These 
were first made available to the official fore- 
casters in the fall of 1959. Although this factor 
may account for part of the increased skill, 
it is also possible that improvement in the 
other tools or differences in the general circula- 
tion may have been responsible. 


It is interesting to note that the official fore- 
casters were able to improve the objective pre- 
dictions when verified in terms of skill score 
and percent exactly correct, but not in terms 
of percent correct within one class (last column 
of table 3). The objective forecasts had an 
over-all score of 82.0 in this category, a slight 
improvement over the official score of 81.4, 
and a decided improvement over the objective 
score of 78.6 during the previous year (table 2). 


Two cases in point are illustrated in figures 
9 and 10. In the first forecast, prepared Nov. 
12, 1959, a large area of below normal on the 
objective (fig. 9 A) was modified to much 
below normal on the official forecast (fig. 9 B), 
in part because of extremely cold air and high 
sea level pressure initially available in the 
polar source region over western Canada, and 
nearly all of this much below normal actually 
verified (fig. 9C). A few days later, on Nov. 
17, the objective prediction (fig. 10 A) called 
or a very large change, from much below to 
much above normal, in the Great Plains, but 
the official forecast (fig. 10 B) correctly modi- 
fied this to above normal in view of the initial 
presence of cold air and snow cover in much 
of the area. A case when the official forecaster 
successfully called for a large change in regime 
despite a poor objective is shown in fig. 11, 
while a case where the official forecaster in- 
correctly modified the objective is illustrated 
in fig. 12. 


OBJECTIVE FORECAST 
NOV, 17-21, 1959 


OFFICIAL FORECAST 
NOV, 17-21, 1959 


MEAN TEMPERATURE 
DEPARTURE OBSERVED 
NOV. 17-21, 1959 


Fig. 12. Objective forecast (A), official forecast (B), and 

observed (C) temperature anomaly classes for the 5-day 

period, November 17—21, 1959. The skill scores were 

22 for the objective and 3 for the official; the percentages 

correct within one class were 92 for the objective and 
73 for the official. 


8. Suggestions for Improvement 


There are several sources of error in the 
objective method as currently applied. In the 
first place, the 36-hour baroclinic forecasts 
(plus the normal) are imperfect representations 
of the desired s-day mean heights for the H, 
period. Even if these forecasts were completely 
accurate for the first 36 hours, there would 
still be errors in estimating H, because of 
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changes in the circulation pattern which ma 
occur during the next few days. A smaller 
error can be attributed to the fact that the 
current 5-day mean surface temperatures (Ty) 
are not perfectly estimated because of inaccura- 
cies in the district forecasts for the next 2% 
days. As a result of these limitations in esti- 
mating values of H, and Ty, the objective 
temperature forecasts tend to be conservative; 
i.e., they are more accurate for the beginning 
than the end of the 5-day period. In fact, a 
separate test was made which showed that the 
objective forecasts explained 40 percent of 
the temperature variance for the s-day period 
centered two days after forecast day (T;), 
“ almost twice as much as the 21 percent they 
explained for the T, period (section 5). 

Even if perfect forecasts of H, and T, were 
available as input for the objective method, it 
could explain at most only an average of 60 
percent of the temperature variance by the 
relatively simple technique used here (section 
4). A more sophisticated approach utilizing 
actual air trajectories and stratifying the data 
might result in greater accuracy. Further im- 
provement will occur when better forecasts 
of the 700-mb circulation for the H, period 
and each of its daily components become 
available. Consideration should also be given 
to the initial thermal distribution over a region 
larger than the contiguous United States, with 
particular reference to abnormalities of tempera- 
ture in air mass source regions. Finally, and 
perhaps most important, the effects of param- 
eters other than 700-mb height and surface 
temperature should be included, such as prog- 
nostic thickness, sea level pressure, snow cov- 
er, ocean temperature, and cloudiness asso- 
ciated with vertical motion. 


9. Conclusion 


An objective method of forecasting 5-day 
mean temperatures over the United States 
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has been developed. For this purpose a statistical 
screening technique was applied to a large 
mass of historical data to produce a fairly 
reliable and stable set of multiple regression 
equations. It has been demonstrated that the 
best numerical forecast of upper level heights 
currently available as input for these equations 
is the 36-hour baroclinic 700-mb prognosis 
regularly prepared by JNWP from 1200 GMT 
data on forecast day. By use of these heights, 
a reasonable level of skill in predicting 5-day 
mean surface temperatures has been main- 
tained over a 2-year test period. 


Since the objective temperature forecasts 
can be prepared in a few minutes with the aid 
of the IBM-704, they should serve as a valuable 
and practical forecast tool. Results during the 
single year that this objective aid has been 
available to the official forecaster are encour- 
aging. It is hoped that continued improvement 
will occur in the future as a result of greater 
accuracy of the numerical height and district 
temperature forecasts, a greater familiarity on 
the part of the forecasters with pitfalls and 
strong points of the objective method in individ- 
ual situations, and refinements in the method 
itself. 
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Abstract 


An earlier theory of steady state, saturated convective currents is modified to include the 
drag of the condensed liquid water on the rising air. The system of equations is numerically 
integrated for several cases in which the entrainment is allowed to assume negative values. 
Results appear to indicate improvement over the earlier model. 


Introduction 


In a recent paper (HALTINER, 1959) a theory 
was presented for a steady state, saturated con- 
vective current undergoing an exchange of heat 
and momentum with its environment by lat- 
eral diffusion, as well as dynamic entrainment 
of environment air. When mass continuity con- 
siderations would result in the ejection of mass 
(detrainment) assuming a constant cross-sec- 
tion area of the updraft, the detrainment was 
taken to be zero; and the cross-section area 
of the current was allowed to increase in order 
to satisfy the equation of continuity. This im- 
plies that when the current begins to decelerate, 
the entire mass slows down and spreads out 
laterally. While this process seems plausible, 
some observations of cumulus clouds indicate 
that detrainment does take place at certain 
times in such convective currents. 


The present paper differs in two respects 
from the preceding theory. Firstly, a specific 
term is included in the momentum equation 
for the drag of the condensed water which is 
carried along by the updraft. Secondly, the 
computations include some examples in which 
detrainment is allowed to take place, assuming 
the mass lost to have zero vertical velocity. 
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List of Symbols 


IK temperature 

w vertical velocity 

q specific humidity 

p Pressure 

0 density of air 

o cross-section area of current 

l specific liquid water content of air 
M = ewo mass rate of flow 

r relative humidity 

D drag of the air on the liquid water 
E external forces per unit mass 

h cloud height 

t time variable 

z height variable 

C, specific heat at constant pressure for air 
c specific heat of water | 

L latent heat of condensation 

Ra gas constant for dry air 

R, gas constant for water vapor 

g gravity force per unit mass 

ye lapse rate of environment 

k, diffusion coefficient 

Subscripts 

e refers to environment 

v refers to virtual 

O initial value 
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The Basic Equations 


Most of the basic equations are similar to 
those of the preceding paper, hence the deriva- 
tions will be omitted except where significant 
change has been made. 

The continuity equation is unchanged and 
is expressible in the form 


dN_ 1 dw 
dz w dz 


1240 80. dl; 
ei odz Rils Ride (1) 


The equation of motion may be written as 
nou Di (2) 


The external forces per unit mass E are the 
pressure and gravity forces, which may be 
combined in the form g(T, - T,.)/T,., and 
the frictional forces. The last of these wi 
include the loss of momentum by lateral 
diffusion and also a specific term represent- 
ing the drag of the liquid water carried along 
by the ascending air. The equation of motion 
for liquid water rising with the air may be 
written as 


dw 
IT -k+D () 


Here Î is the mass of condensed water per 
unit mass of air, and D represents the drag of 
the air on the / grams of liquid water. Thus 
— D, which may be obtained from Eq. (3) is 
the retarding drag of the liquid water on the 
air. The loss of momentum by the convective 
current per unit mass per unit time through 
lateral diffusion to the environment is represent- 
ed by a term of the form — k, (w - w.), where 
k, has dimensions of sec~!. Thus the combined 
“frictional” forces per unit mass of air are 


uw) Ile) (a) 


It follows from Eqs. (2) and (4) that the equa- 
tion of motion for the rising air is 
dw = 1 x, 1 
co 


—(w — we) m — ky (w-w.) - 
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Utilizing the assumption of steady state leads 
to the form 


dw g/[T,- 
(x +2) sf 


-1) (ww) by 
(5) 


The thermodynamic equation is essentially 
unchanged from the previous theory, namely 


d dT dN 
= de [c,(T - Te) +L (q - q)l ES 
k 
_ i E (T- T.) +L (q - 4) | = 
dT ” 
ge (6) 


The second term on the left, which represents 
the heat energy given up by the liquid water 
as it ascends with the air, is only a few percent 
of the first term on the right, namely, the 
enthalpy change of the air. The former is in- 
cluded here mainly for consistency since the 
effect of the drag of this liquid on the air is 
included in the momentum equation. However 
since it was considered negligible and omitted 
in the previous results, it will also be omitted 
in the computations for the present model in 
order to afford a more accurate measure of the 
more important effects of the drag term and 
detrainment. It might also be pointed out here 
that in nature the water droplets acquire a 
downward velocity relative to the air and 
Eq. (3) would not apply in general. 

The Clapeyron equation for the saturation 
vapor pressure together with the definition 
of specific humidity gives 


Id E dT | 8 
ade Roi" de or. (7) 


Since detrainment is being considered, the 
law for the liquid water content / must be 
piecewise defined as follows : 


(i) during entrainment, 


dl d d k 
Re ge ec 
(ii) during detrainment, 
dl da k 
en ge Ga) (9) 
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Using the definition of virtual temperature With appropriate initial and environmental 


T, = T (1 + .61 q) (10) 


it may be shown that 


dT, 


v dT 
Je aE + 12.5 DT (11) 


In this study the cross-section area of the 
current will be assumed constant. Thus with 


(12) 


and some minor approximations, the preceding 
equations may be expressed as the following 
system which is in a convenient form for the 
numerical integration. 


ar fz 4 
a fal gern] x 


esl eV 2 ] 
w ee HR 


ky Cp gly &4 
+l q- qt 27-7) + ~ 


do =o 


ET Ryle 
Lq _s_d (G+r254, 
RME dE ER 


- te eg 
La Ye er n]) 
dN If g (= le ) ky g 


dz 2Lw? Le 


SIT, 
ar: Z| (14) 
dw I dN g 
Dep -(v Pre pe 
D Je 
(7 T -1)] (15) 
dg Lg ates eq Z 
DR Pee RTs (36) 
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conditions this system of equations may be in- 
tegrated by numerical methods to yield val- 
ues of T, N, w, q and |. The factor (1 +)! 
in Eq. (15) may be replaced by 1 with a high 
degree of accuracy. With this minor approxi- 
mation, the above system differs from the 
equivalent system in the previous model H [2] 
in two respects. Firstly, because of the dra 
of liquid water, the term [(T, - T,,)/T,.| - | 
appears in several equations above replacing 
(T, — Tye)/ Tye in the previous theory. Secondly, 
because some of the computations to follow 
include detrainment, Eq. (17 b) must be used 
to determine the liquid water content after 
the detrainment begins. In the computations 
to follow, it will also be assumed that w, = 0, 
as in the previous computations. 


Results and Conclusions 


The system of differential equations 13—17 
was integrated by the Runge-Kutta-Gil meth- 
od (GILL, 1951) on one electronic digital com- 
puter. Four sets of initial and environmental 
conditions were used, each identical to a set 
of conditions already computed for the previ- 
ous model. 

For purposes of comparison, cases (a) and 
(b) have the same initial and environmental 
conditions; however the term representing 
the drag of the liquid water - D was omitted 
in case (a), while retaining the detrainment 
feature of the present model. Cases (a), (c), 
and (d) include both water drag and detrain- 
ment. The conditions corresponding to the 
various cases are as follows : 


Case (a) 
D COM CN TE FIIR 
Wy = Imjsec, w, = 0 


Ve 7 CJI00N 1 = 80% 

ky = 001 sec: +, D> = 0 
Case (b) 

same as case (a) except D + 0 
Case (c) 

same as case (b) except that y. = .6° C/100 m 
Case (d) 


same as case (b) except that r = 90 % 


Figures 1, 2, 3, and 4 present a comparsion of, 
respectively, temperature excess over environ- 
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Figure 1. Temperature excess over environment as a 
function of height for Case (a), Case (b) and after HALTINER 
(1959). 


ment AT = T-Te, fractional increase of 
rate of mass flow M/M,, vertical velocity w, 
liquid water content and also the corresponding 
results of the previous theory (H). Cases (a) 
and (H) should be the same up to the point 
where detrainment begins, however a small 
systematic difference may be noted. This minor 
difference is due to slightly different computa- 
tional procedures, primarily in the determina- 
tion of the environmental parameters. 

Figure 1 shows that the differences in tem- 
perature among the three models generally 
are only of the order of a few tenths of a 
degree, which would make it rather difficult 
to indicate the preferable model on the basis 
of temperature observations in clouds. Figures 
2 and 3 show that permitting detrainment, 
while omitting the drag of the liquid water 
(case (a)), gives the greatest cloud height and 
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Figure 2. Fractional mass change M/M, a as function 
of height for Case (a),Case (b) and after HALTINER (1959). 


vertical velocity. This is to be expected since 
the ejection of mass (detrainment) tends to 
accelerate the remaining mass of the ascending 
current. On the other hand, the inclusion of 
the additional frictional drag of the liquid wa- 
ter tends to reduce the cloud height, vertical 
velocity and, in turn, the total cloud mass, the 
latter being represented by the term M/M,. 

It is interesting to note that the apposing 
effects of detrainment and the retarding drag 
of the liquid water in case (b) result in nearly 
the same cloud height as (H) where both are 
omitted. However, the maximum vertical 
velocity attained in case (b) is significantly 
less than that of the previous theory (H), or 
case (a) which does not include the drag of the 
liquid water. This is to be expected since re- 
tarding drag of the liquid water begins im- 
mediately, whereas the accelerating effect of 


Tellus XII (1960), 4 


SOME FURTHER RESULTS ON CONVECTIVE CURRENTS 


90 


80 


70 


o 
(e} 


a 
O 


HEIGHT (HUNDREDS OF METERS) 
oi > 
re) O 


20 


w (METERS/SEC) 


Figure 3. Vertical velocity as a function of height for 
Case (a), Case (b) and after HALTINER (1959). 


the detrainment does not begin until after the 
maximum vertical velocity has been attained. 

Differences in the liquid water content of the 
cloud are rather small, averaging about 10 %. 
Table 1 summarizes the results of the com- 
putations. 

In the previous paper computed values were 
compared to some observational data from an 
actual cloud. In general the agreement between 
observed and computed parameters was fairly 
good; however the computed vertical velocity 
was somewhat in excess of observed values. To 
the extent that the modifications reduce vertical 
velocity, while comparing equally favorably 

in other respects, the modified model appears 
to be an improvement over the previous theory. 
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Figure 4. Liquid water content as a function of height 
for Case (a), Case (b) and after HALTINER (1959). 


For a more complete understanding of 
convective phenomena further work is neces- 
sary in regard to the cross-section area of the 
convective currents, the environment vertical 
velocity and the horizontal distribution of 
convective cells. These questions will require 
the introduction of horizontal variations and in 
essence a three-dimensional analysis of convec- 
tion. Morcover, more knowledge must be 
gained about the initiation of the downdraft 
in a large convective cloud, as well as about 
the initial development of the convective 
current itself. In connection with the latter, 
some computations by Markus and Wirt 
(1959) show a bubble-like structure during 
the development of a convective element. 
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Table 1. Comparison of maximum values of AT, M/M, I, w and cloud 
height, h, with maximum values after Haltiner (2). 


SE m nn u 


Case Parameter (H) Case (a) Difference 


M/M, 8.13 8.56 43 
w (m/s) 14.12 15.52 1.40 
1 (g/kg) 6.0 6.4 i 
h (m) 7900 8750 850 
(H) Case (b) 
AT 1.26 1.60 34 
M/M, 8.13 6.22 1.81 
w 14.12 10.13 3-99 
1 6.0 6.83 83 
h 7900 755° 350 
Case (a) Case (b) 
ATP 1.39 1.60 Pai 
M/M, 8.56 6.22 2.34 
w 15.52 10.13 5.39 
ı 6.4 6.83 43 
h 8750 7550 1200 
(H) Case (c) 
AT — .0I .97 .98 
MM, 2.88 1.92 .96 
w 4.07 2 1.95 
l 5-9 4-2 1.7 
h 4759 4500 250 
(H) Case (d) 
ANT 2.47 2.76 29 
MM, 11.80 9-75 2.05 
w 22.28 17.78 4.50 
1 Tak 77 .6 
h 9500 9250 250 
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Abstract 


A liquid of mean depth H is rotating at constant 
angular velocity (@) inside a cylindrical annulus whose 
radial width is large compared to H and small compared 
to the mean radius. An equilibrium Ekman flow is 
then produced inside the boundaries by (for example) 
symmetrically forcing liquid into the annulus at the 
outer rim and removing the same quantity at the inner 
circumference. In response to the radial pressure gradient 
a geostrophic azimuthal velocity ([J) is established above 
the inflowing viscous boundary layer; the latter extending 
to a characteristic height hp = (v/w)'/2 above the bottom 


surface of the annulus. This is a study of the shearing 
instability of an Ekman flow whose components are 
uniform in both horizontal directions. 


We first discuss the anticipated similarities and dif- 
ferences of the stability of this “rotational” boundary 
layer and their counterparts in hydrodynamical flows 
which have no Coriolis force. Thus we note that the 
available kinetic energy in the boundary layer of the mean 
Ekman flow can be released by ‘“‘pure-boundary” modes 
whose characteristic dimensions are of order hy, and 


whose onset condition only depends on a boundary layer 
Reynoids number R = Uhzv”!. This paper establishes 


the theoretical possibility of a qualitatively different 
energy releasing process, apparently confined to gyro- 
scopic boundary layers, by “body-boundary’’ modes and 
suggests that for certain values of the Taylor number 
(E =wHy-1) these may occur at lower Reynolds 
numbers than the pure boundary species. The energetic 
picture of this instability is as follows: Pressure fluctua- 
tions in the mean boundary layer, whose radial wave- 
length is much larger than h,, but much smaller than H 
and whose azimuthal wavelength is much larger than H, 
are partially (quasi-hydrostatically) transmitted vertically 
where they are accompanied by quasi-geostrophic azi- 
muthal components of perturbation velocity and ageo- 
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strophic radial components. The coupling or this interior 
component with the boundary component of the per- 
turbation influences the energy which is transferred to 
the latter from the radial (i.e. cross-isobaric) component 
of the mean field, and more than compensates for the 
burden of driving the perturbation against the stabilizing 
viscous forces in the interior of the fluid. Perhaps the most 
noteworthy aspect of these modes is that they can trans- 
port mean energy and momentum from the deep interior 
into the boundary, thereby making it available for 
further dissipation or to replenish the mean radial energy. 
The pure boundary modes, on the other hand, can only 
alter the geostrophic momentum in the interior by pro- 
ducing horizontal convergences of the mean flow in the 
friction layer. 

We have used two slowly “converging’ asymptotic 
expansions to obtain solutions of the characteristic value 
problem and have supplemented these resulis with error 
estimates to indicate the quantitative limitations in the 
range of physical interest. For large values of the Taylor 
number a disturbance of fixed horizontal wavelength 
will be unstable for values of the Reynolds number 
occurring in discrete ranges which are separated by other 
ranges in which the disturbance transfers energy to the 
mean flow. These unstable ranges or “scales’’ correspond 
to the set of eigenfunctions having different vertical struc- 
ture. On the other hand for fixed supercritical Reynolds 
number there is a preferred radial wavelength in the nth 
scale which is given by 2nm, Eh, where m, ap- 
proaches a numerical constant as E> oo. The correspond- 
ing minimum Reynolds number for the marginal stability 
of the nth eigenfunction is R, = x El, where x, is a 
numerical constant which only depends on the integer 
n, in the asymptotic limit E> co. For the “largest scale 
of motion” n= I, this asymptotic limit is easily reached 
under laboratory conditions (e.g. E = 2,500). Because 
of the limitations of the second asymptotic expansion we 
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can only estimate the numerical constants (#3, m1) as: 
m, = 10°, #,<11. Thus, for example, when E= 2,500 
the Ekman flow should not be stable for Reynolds 
numbers as large as 8 x10’. For larger values of E the 
Reynolds number which is necessary to destabilize the 
n= 1 scale increases. However smaller scales of motion 
(larger n) are preferred as E increases since it is found that 


Anti Xp in the asymptotic limit E— oo. The present 


theory is unable to determine that finite value of E 
when the second scale of motion becomes preferred 
ver. 

It is believed that a manifestation of this instability may 
be produced in the laboratory by causing a small decrease 
in the rate of rotation of a cylindrical vessel containing 
a shallow layer of water. Under certain conditions the 
decaying relative motion has the form of a time-modulated 
Ekman flow. By means of a permanganate solution placed 
in the boundary layer, the breakdown of this symmetric 
flow and the formation of thin rolls has been observed. 
We have also correlated experimental measures of the 
average rate of decay of relative angular momentum 
with a detailed laminar calculation (not included herein). 
For the converging Ekman flow, starting at a value of 
E of about 2,500, the measurements indicate a syste- 
matically greater mean decay rate than one would expect 
if the motion were laminar and symmetric. For the range 
of Reynolds numbers which were used in the experiment 
this critical point varied substantially less with changes 
in Y (or R) than with changes in (or E). From these 
qualitative and quantitative observations we have con- 
cluded that the energy source for the instability is in the 
mean Ekman flow but that it cannot be entirely ex- 
plained as a boundary layer phenomenon. It is suggested 
that the body-boundary modes discussed in the instability 
theory provide the mechanism for the increased transfer 
of angular momentum from the interior of the fluid 
to the lower boundary. 


1. Introduction 


In 1905 Exman (described in Lams (1945) 
and numerous meteorological texts) published 
his theory of the equilibrium current which 
results from a balance of viscous, Coriolis, and 

ressure gradient forces in an attempt to 
describe the motion of a homogeneous ocean 
as the result of the stress exerted by the atmos- 
pheric winds. Although the insubstantial nature 
of the Austausch coefficient need no longer be 
dwelt upon, the physical process which the 
theory attempts to describe is undoubtedly of 
major importance for the transport process in 
the surface layers of both planetary fluids and 
to the general circulation problems to which 
they are ultimately coupled. Such flows as are 
predicted by Ekman’s theory may, hôwever, 
be realized in the laboratory by reducing the 
Reynolds number until laminar motion is 
obtained. Granting the relevance of laboratory 
models to the geophysical phenomenon we 
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must pause to assay the utility of the problem 
of laminar instability, with which this paper 
deals, for the study of fully developed tur- 
bulence. Despite the gap in the non-linear 
scale which separates these two, several hydro- 
dynamicists have worked in the belief that 
they are intimately coupled at least in those 
aspects pertaining to the transport problem. 
A recent theory of the mean velocity profile in 
turbulent pipe flow by Markus (1956) aban- 
dons detailed energy.and momentum balance 
on the grounds that there are a multitude of 
solutions to the Navier-Stokes equations. The 
realized or most stable solution is conceived of 
as a discrete set of independent modes which 
draw their energy from, interact through, and 
consequently determine the mean field. This 
generalized eigenvalue problem and the quanti- 
fication of the idea of “most stable solution’ 
leads to a theory of the momentum transport. 

The reader who is aware of the historical 
development of the theory of parallel flow 
instability will not expect the present paper to 
bring the Ekman problem to this state of 
development. The aim has been to explore 
the nature of the energy releasing process and 
to demonstrate by an asymptotic theory the 
possibility of a mode which is apparently 
confined to gyroscopic boundary layers. The 
mechanism which is discussed would appear 
to be of geophysical interest for the following 
reason. One usually conceives of the turbulence 
which originates at the boundary of the carth 
as effecting the flow in the interior of a 
statically stable atmosphere or ocean by means 
of large (synoptic) scale convergences in the 
boundary layer flow. This study establishes the 
physical basis for an additional and more 
direct type of coupling. 


2. Statement of the problem 


Fig. 1 is a schematic picture of a conceptually 
simple Ekman flow, and one which may be 
approximated in the laboratory. The section of 
water which is shown, is at some large distance 
from the axis of rotation of the smooth and 
rigid bottom surface. The homogeneous fluid 
with a free surface at height H is maintained 
in a horizontally uniform and steady flow 
relative to the bottom by means of the radial 
(y‘) mass transport. Either this quantity or the 
geostrophic azimuthal current (U), which is 
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Figure 1. Schematic diagram of the mean Ekman flow 
and the body-boundary perturbation. The section is 
shown at a large radial distance from the axis of rotation 
(@). Vo is the “radial’’ (y’) or cross isobar component 
of the undisturbed current. (/ is the amplitude of the 
azimuthal (x’) component at or near the free surface. 
The radial half-wavelength of the perturbation is 
m%HE~/|sm-1 where m is an order unity wave number, 
H is the height of the free surface and E is the Taylor 
number. ve denotes the radial disturbance velocity in 
the interior of the fluid and the vertical concentration of 
the streamlines at the rigid bottom illustrates the boundary 
component of the mode. 


produced by the hydrostatic transmission of 
the radial pressure gradient to the free surface, 
may be regarded as a primary known quantity 
in the analysis of the stability of this laminar 
Ekman flow. 

Since the mean field is independent of the 
horizontal coordinates the non-linear terms in 
the equations of motion vanish and we can 
obtain an exact solution of the equations of 
motion involving a balance of viscous, Coriolis, 
and pressure gradient forces (e.g. Lamb sec. 
334 a). This solution for the mean field is only 
a function of the vertical coordinate 2’, and 
shows that the ageostrophic motion is con- 
centrated in a viscous boundary layer whose 
thickness is of order hg= (v/w)': where « is 
the rotation rate and » is the kinematic mo- 
lecular viscosity. Dimensional considerations of 
the Navier-Stokes equations with hg as the unit 
of length, =! as the unit of time and U as the 
unit of velocity lead to the conclusion that the 
instability of this flow can only? depend on 
a boundary layer Reynolds number (R= 


1 In accord with the neglect of the horizontal variation 
of the basic current we neglect the variation in height 
of the free surface. More precisely we have assumed a 
rotational Froude (F= U?/gH) number which is small 
compared to one. In this attempt to isolate the mechanism 
under discussion, we have also neglected the rim boundary 
currents as a source of instability. 
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=U(vr/o)'®v!) and an Ekman or Taylor 
number (E= @H?y-1) which measures the 
ratio of total depth of the fluid to the thickness 
of the boundary layer. 

Since the basic field gradients are confined 
to the boundary layer, the source of energy 
for a growing disturbance, or at least one of 
infinitesimal amplitude, must be in this region. 
There is then a tendency to classify the stability 
problem with other boundary layer shear flows. 
In the flow over a flat plate (Blasius), for 
example, the experiments of SCHUBAUER and 
SKRAMSTAD (1947) indicate that the total depth 
is not a significant parameter in the instability 
as long as it is much greater than the laminar 
boundary layer thickness. The measured value 
of the critical boundary layer Reynolds 
number is above 400. We also cite the case of 
the stability of the three-dimensional boundary 
layer produced by a circular disc rotating with 
angular speed w in an ambient atmosphere of 
viscosity v. The boundary layer thickness for 
this centrifugal fan is also (y/w)’!. GREGORY, 
STUART and WALKER (1955, p. 190) cite a 
range of 550—680 for experimentally deter- 
mined value of Ra. In both examples the 
preferred horizontal dimension is proportional 
to the boundary layer thickness. If the unstable 
modes for the Ekman profile are also confined 
to the boundary layer, then the criterion on- 
ly depends on R but not E. The analytical 
problem is more complicated than other shear 
flow instabilities because the mean field is not 
a parallel flow, and in addition it has a unique 
type of structure which is related to the Coriolis 
force. Nevertheless, it is a reasonable presump- 
tion that the Ekman flow will be unstable to 
boundary-like perturbations at Reynolds num- 
bers of the same order as in the case of the 
centrifugal fan. 

Preliminary experiments, to be referred to 
in Section VIII, suggested that the Taylor 
number was also a significant parameter. If 
this be the case, then one must admit a class of 
modes which extend throughout the depth of 
the homogeneous fluid and which satisfy the 
boundary conditions at the upper surface. Con- 
sidering that no kinetic energy is available to 
the perturbation in this interior region, we 
may well ask if it is not more economical for 
the disturbance to confine itself to the bounda- 
ry. The following asymptotic theory suggests 


two partial answers to this question. One is 
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that the coupling of the so-called interior and 
boundary components of the perturbation 
increases the transfer of energy to the latter 
from the boundary layer of the mean flow, and 
more than compensates for the burden of 
driving the interior component. Secondly, 
disturbances on this scale perform a function 
which the pure boundary modes cannot. They 
are able (in the finite amplitude stage) to 
transfer mean energy from the remote interior 
and make it available in the Ekman layer for 
further dissipation. At large Reynolds numbers 
both species may be present. 


3- Qualitative Remarks on the Energetics 


In order to motivate the unusual dimen- 
sionalization of the equations of motion and 
the subsequent asymptotic expansion in the 
Taylornumber, we make the following prelim- 
inary remarks concerning the free dissipation 
rate of disturbances in a rotating fluid. By this 
we mean that the energy source which drives 
the mean relative motion in Fig. 1 is removed 
and the fluid is in solid rotation except for a 
small initial perturbation V’, which decays in 
time (f’) due to the action of molecular viscosity. 
A class of disturbances which retain their shape 
during this decay process (i.e. the charac- 
teristic functions) may be found by writing V’ 
as exp| —A’t’ +2ziy'/L,] times a function of z’, 
where Ly is the wavelength of the disturbance 
in the “radial” direction and A’ is the decay 
rate. Since the discussion has been simplified 
by only considering those modes which are 
independent of x’, the decay rate will only be 
a function of (w, », H, L,). We assert that 
under certain conditions (to be determined) 
this eigenfunction has a boundary layer struc- 
ture, which in fact is similar to the classical 
Ekman profile for a mean flow that is both 
steady and independent of the horizontal 
coordinates. It consists of interior components 
(ue, Ve) which are independent of 2’. The 
azimuthal velocity u. is in quasi-geostrophic 
balance while the radial velocity v, is associated 
ue 
Os 
Since the interior components of the mode 
cannot satisfy the no-slip condition at z’= 0 
we must add a viscous boundary layer compo- 
nent (#, v5) extending to a charactesritic height 
hz= (vo) above the bottom surface. A 


with interior momentum changes (20)! 
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vertical section of-the total perturbation is 
schematically indicated in Fig. 1. Almost all 
of the perturbation kinetic energy is in the 
interior, while almost all of the dissipation is 
in the boundary layer. The latter region is 
almost in equilibrium during the decay process 
as the result of its pressure coupling with the 
interior component, which in turn decreases 
the local kinetic energy by inviscid (almost) 
radial displacements. The perturbation kinetic 
energy averaged over one wavelength is then 
of order #2H. Since the boundary layer 
amplitude is of the same order as u, it follows 
that the viscous dissipation due to vertical 
shear is of order v(hz2)uhr. On the other 
hand, the rate of dissipation due to lateral 
shear, a body effect, is of order »L,?wH. 
These two dissipation rates are of the same 
order when Ly, ~ (Hhg)’*= E'«hz= EH. 
When L, is of larger order of magnitude than 
hgE!* the interior viscous force is negligible. 
Since total dissipation must equal the rate of 
decrease of kinetic energy 4’u2H = vurhz' + 
+vL 42H. We therefore conclude that when 
L, is of order HE (or greater), 2’/w~ 
= vo hg H-! = Eh, Thus the horizontal dimen- 
sion of the mode may be much less than the vertical 
dimension and still have no influence on the 
dissipation rate, a situation which is in marked 
contrast with the analogous problem in a non- 
rotating fluid of depth H where the decay rate 
increases as L,? and is independent of H 
for L, > H. (For a more rigorous justification 
of these conclusions see the remarks following 
eq. (21) in Sec. IV.) 

Let us now return to the stability problem 
where the question is one of balance of transfer 
and dissipation of energy. Although the mean 
field will affect the shape of the characteristic 
function, perhaps the previously discussed 
relation between dissipation and horizontal 
wave-length will be qualitatively unaltered, 
i.e., the dissipation increases slowly as the radial 
length of the mode decreases. At the same time 
the vertical component of the perturbation 
velocity increases relative to the horizontal 
component, for reasons of mass continuity. We 
therefore expect a relative increase in the 
Reynolds stress, as the horizontal dimension 
decreases, and consequently an increase in the 
transfer of energy from (or to) the mean flow. 
The qualitative idea that the transfer increases 
more rapidly than dissipation and eventually 
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overtakes it at some wavelength of order 
HE: is explored in a more formal fashion 
below. 


4. Asymptotic expansion in the Taylor 
number 


We therefore choose a unit of time (w~1E’/) 
equal to the free viscous decay time of these 
modes whose horizontal dimensions are much 
larger than the Ekman depth. With a unit of 
length equal to H and a unit of velocity (U) 
equal to the magnitude of the basic geostrophic 
flow, the equations of motion for the growth 
of infinitesimal perturbations (V’) are: 


ON 
ae _E1v?V'+28 x V'+ 


Sab (Vo ev PV eV) = = vp (1) 


v-V =o 


Er REV 2) 
San 3 


where x is to be regarded as a new independent 
parameter in lieu of the Reynolds number, 
and the scaling parameter f is to be determined 
in such a manner that it leads to a consistent 
asymptotic expansion at large E. Z is a unit 
vector along the axis of rotation, and vp’ 
is the normalized pressure force. The Ekman 
depth is now E-‘, and the basic field may be 
denoted by V,= V,(zE":), having a compo- 
nent V, in the y direction which decreases to 
zero at large values of zE': and having a 
component I+U,(zE':) which approaches 
unity. 

Since these fields are independent of the 
horizontal coordinates the eigenfunctions of 
(1) will be sinusoidal functions of (x, y). Let 
us now introduce the following scaling trans- 
formations for the (x, y, z) components of 
the velocity V’= (u’, v’, w’): 


fu, vw’, p} = exp (At + iymE™ + ixlE-*) x 


WU, Vows. DE) 


The reason for the choice of the y scaling in 
the above equation should be apparent from 
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the remarks in Sec. 2, and we are free to assign 
values to a, b, c, f which will lead to a 
simplification of the eigenvalue problem at 
large E. Since the amplitude in a linear per- 
turbation problem is arbitrary we may take 
the amplitude of the interior azimuthal velocity 
component (u,) as unity. If we assert that the 
y component of the pressure gradient is of the 
same order, then since m is defined to be inde- 
pendent of E, c=!/,. For the balance of 
azimuthal (x) interior forces we have indicated 
in Sec. 3 the importance of the Coriolis term 
proportional to v’ and the inertial term pro- 


2 


portional to &-1. Since the latter is of 


ou 
ae 
order E~‘ in the interior then the amplitude 
of v’ should be chosen of order E-‘ in the 
interior. The azimuthal scaling factor, a= %, 
is then determined by requiring that the pres- 
sure gradient in the x direction is at most of 
the same order as the interior Coriolis force. 
In order to have a coupling between the 
boundary layer and the interior Aw’/dz in 
the latter region must be of the same order 
as the horizontal convergence, and this leads 
tob= 4. 

The sole remaining scaling parameter is f 
and this will be chosen to make the transfer 
forces (like E~fw'IV,/dz) of the same order 
as the viscous force (E-1 v ?V’) in the boundary 
layer. With our previous choices the former is 
of order EEE"; while the latter is of 
order E-1E+1, and therefore we take f= 4. 
In short, we introduce the transformation: 


{u’, Ti w’, ps = exp (At ar iymE'!: + 
+ ixlE Ts) {u, v, E“hw, Ep} 


R=E ix (3) 


where 2, u, v, w, p are functions of the new and 
independent parameters (x, E,/,m) and z. 
Equations (1) and (2) then become: 


i a? à : 
+24 + 2] im Vov + Ehilo(x + Us) + (4) 


found : 
+ E- wo Vo(zE | = — imp 
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(Dy A ES SEA (x — mE: — PE) u = 
dz? 
—2v+% | imVou +E~hil(1+ U,u+ (4) 


+ E y U | = -ilpEr 


2 
() Eh (= MER PEI + 


dz? 
+ x[imV. w+ E~hil (1 +U )w] = — 2p 
> 5 Oz 
en = - Elsimv - ilu (5) 


Equations (4) are of course equivalent to the 
original equations (1) and (2) since we have 
merely re-defined variables. Now however we 
introduce the assumption that for arbitrary 
but fixed values of (x, J, m) the solution may 
be expanded in the small number Eh as 
follows: 


u =teo(H, 1, m) + En. (2,%,1,m) + 
+ Elu, [(1 -z)ER] +... (6) 


+ Uno (C, x, 1,m) + E hum +.... 


pe ES hy, (x1, m) +B We (2, elm) + ove: 


+ Vo9(C, x, 1, m) + Eh vy (C, x, L m) +... 


(7) 
A= (x, m) + EA, (2, lm)... (8) 
C=zE 


The subscript e denotes the body component 
of the perturbation, which describes the 
horizontal velocities in the region of vanishing 
mean field gradients. The leading term in (6) 
and (7) are independent of 2; to in fact 
being equal to 1 while the leading term for the 
radial interior velocity is of order Eh, for 
reasons mentioned above. The boundary com- 
ponents of the mode, denoted by subscript b, 
are rapidly varying functions of z and which 
like U, and Vy decrease to zero above the 
Ekman boundary layer. Their leading terms 
are of order unity since viscous forces are 
of the same order as inertial forces in the 
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boundary layers. The last term (uz) in (6) 
identifies a boundary component at the free 
surface which is necessary, in high order 
iterations, to nullify the surface stress produced 
by the component #1. Substituting (6) and (7) 
in the continuity equation (5) and integrating 
gives the following expansion for the vertical 
velocity as a function of z or {=zE': 


w= -iz (mv + Ingo) — jEh f [mva + In.]dz 
. 0 
5 iq 
—imf vyo() dé —iE "hf (mvi + luyo) dl (9) 
0 0 


If the pressure perturbation at the undisturbed 
surface (z=1) is denoted by po then integra- 
tion of the third equation of motion gives the 
following relation between pressure and veloc- 
ity components: 


p = Peo 3 Ehpa(z) 725 Mage 
+ E~ py (£) Sic .... 


(10) 


where 


palz)= - “(dg +m? + ilx) (mves + lues) (I — 2?) 
(11) 


por(C) = — umvyo(C) + imvy9(E = BE’) 


El: 


+ 4m [ Voll) de’ fy(L’) de” (12) 
¢ 0 


The velocity functions (we, ve) and (up, vs) 
are now determined by requiring that the 
former satisfy the dynamical equations (4a) 
and (4b) in the interior of the fluid (2> Eh) 
where V, and U,>0, and by requiring that the 
latter satisfy the boundary layer equations 
formed by substracting the interior compo- 
nents from equations (4a) and (4b). When the 
coefficients of the terms containing the smallest 
power of E~‘/ are set equal to zero, we get 
the following equations corresponding to (4a) 
and (4b) for the zero order interior compo- 
nents : 


(a) 2e + iMPeo = O 


(b) Ap theo + M? ten — Wen + ul + ilpeo = © (13) 
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One then finds the zero-order boundary layer 
equations by substituting the expansions 
(6)—(12) in (4a)—(4b), cancelling terms cor- 
responding to (13), and stretching the vertical 
coordinate by the transformation z=EE-"k, 


— v6o(C) + 2460 + imx [Vo (2) vb0 - 
[a 
= VoS Vi de] = O 
0 
= yo (2) — 2060 + ime [Vo (no + eo) — 


= Ups vio de] = © 
0 


(14) 


The sum of the two velocity components 
must vanish at € =o and therefore the follow- 
ing zero order boundary condition must be 
satisfied: 


ty (0) = — ten 


(15) 


bo (o) = © 


The two components of the perturbation are 
also coupled by the boundary condition on w 
at the upper surface. Thus when (9) is evaluated 
at z=ı (or €=E*). and the result is set equal 
to zero we get (again to zero order): 


(16) 


oo 
MVeq + leo = — mf Vode 
0 


Equations (13)—(16) are sufficient to specify 
the zero order eigenfunctions and growth rate 
once the amplitude of the perturbation is given. 
Since this quantity is arbitrary, we take u,9 = I 
for convenience and refer all perturbation 
amplitudes to this unit. Then when the pres- 
sure term is eliminated from equations (13) 
and v,, is eliminated by (16) one gets the 
following expression for A, in terms of the 
integrated boundary layer velocities: 


Ag = —m? -ile - ns (C)de 


0 


Red, = - m? + 2 Red (co, mx) 


5 
where A(C, mu) = — Svwld)dE 
0 
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The new dependent variable & is the vertical 
velocity due to horizontal convergence in the 
boundary layer and a single equation for it is 
obtained by combining the two equations in 
(14) as follows: 


29 = - D + imu (Vod' - Vo'¢) 
— Ugo +26 + imx [Vo (ug + 1) + Usd] =0 
(20) 


or 
BY —imx|2V (6) 6" + Vob" - Vo'd'| +46 + 


+ mb (EU, +Vo') - mV (VD — Vod) = 


(21) 


= —imxV, 


The asymptotic expansion in the Taylor 
number has reduced the characteristic value 
problem to the solution of equation (21) which 
contains only a single parameter mx. The 
important coupling with the body mode lies in 
the inhomogeneous boundary condition. The 
stability characteristics are obtained by evaluat- 
ing this solution (6) at {= and substituting 
in (18). It is noted that when mx=o (no 
mean field) eqs. (14) have exactly the same 
form as the non-dimensional equations for an 
undisturbed Ekman boundary flow. Their 
integration is then readily accomplished and 
one finds d(®, 0) = —1/, which corresponds 
to a decay factor exp —{(1+m?). For m<ı 
these “Ekman-like”” modes have a decay rate 
independent of wave number, which agrees 
with the dimensional argument in Sec. 3. 

It is a consequence of eq. (18) that if 
the real part of (00, mx) is greater than zero 
then the flow must be unstable to infinitesimal 
disturbances. for sufficiently large Taylor num- 
bers and sufficiently large Reynolds numbers. 
The marginal stability curve is obtained by 
setting (18) equal to zero and solving for x as 
a function of mx, which then results in: 


(mx)? 
_ 2Red(co, mx) 


(22) 


lim: x? 
E->oo 
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Notice that the wave number / does not appear 
in (22). This degeneracy means that as long as 
the azimuthal wavelength is much larger than 
the depth of the fluid there is no preference 
amongst all such infinitesimal disturbances. 

Eq. (18) may also be interpreted as an energy 
balance for the interior motion, where the 
first term on the right is the viscous dis- 
sipation and the last term represents the energy 
which is transferred at the top of the boundary 
layer by the action of the perturbation pressure 
force. We have also obtained an expression 
for the first correction to this growth rate by 
setting the coefficients of the next largest 
power of Eh equal to zero. This is listed 
below for the purpose of obtaining the limits 
of validity of the foregoing expansion. Thus, 
correct to order E-1 


A= -m? - ilk +26(00, mx) 


oo 


I co 
+ eal me = [ vd = 2 [vu at | 
0 


0 


From the first term in the brackets we estimate 
that the zero order dissipation (m?) is in error 


by a factor | 1 2 Eh Reds , and since 
2 3 


2 
m ar 
Rep = Bis marginal stability the error factor 


is of order [ ai = Bava me |. Form < req. (18) 


or (22) is quite adequate for any Taylor 
number of practical interest. For m>1 (but 
independent of E) the Taylor numbers which 
are required by this expansion soon exceed 
laboratory possibilities. 


5. Solution of the characteristic Equation (21) 
for |mx|>1 


The analytical complexity of (21) is such 
that we must again resort to asymptotic 
techniques. We could expand the solution 
as a power series in small values of mx or, 
alternatively, in large values of |mx|. We 
choose the latter technique in order to demon- 
strate conclusively that & (co, mx) is positive 
for some values of the parameter. We believe, 
however, that the vicinity |mx| ~1 which 
is beyond the reach of the following ex- 
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pansion, is of most interest for computing 
minimum critical R under laboratory condi- 
tions. 

The mean Ekman boundary fields (Up, Vo) 
are determined by the hydrostatic balance of 
viscous, Coriolis, and pressure gradient forces. 
Since & is measured in terms of hg the appro- 
priate equations are: 


d yur 
gel + Vo)=0 
re U, =o 


(23) 


0 =V,(0) =1 + U,(o) = (co) = U,() 
Using the first of these equations to eliminate 
U, in (21) we get the fifth order inhomogene- 
ous equation 


$Y - imu (aVib” + Vid - Vib’) +46 
- mx2V_(Vod' —Vob) = —imxVo(C) 
(24) 


The additional relations in (23) also lead to a 
unique solution for V,, viz. Vy=e- sin£. In 
adapting the asymptotic techniques that have 
been evolved for the Orr-Sommerfield equa- 
tion (LIN, 1945) to the present problem, we 
will see that artificial singularities are intro- 
duced wherever V,(£) =o. Since successive 
maxima and minima in the realized profile are 
in the ratio of the small number e“, let us 
first consider the formal problem of evaluating 
the integral p(c, mx) for an arbitrary func- 
tion V4(¢) which is positive for all values of & 
except at €=o and Ë =, where V,(o) = 
= V,(~) =o, V;(0) = 1. It is anticipated that the 
desired integral property may be approximated 
by fitting some such positive definite Vy to the 
realized profile e~ sin € and that the associated 
error may be estimated from the r.m.s. de- 
parture of the two profiles. Although such 
smoothing may eliminate some modes, it is 
unlikely that the small ripple component could 
completely stabilize a mode which is found 
to be unstable by this approximation. In 
particular we shall consider the function: 


Vo(C) =e sine, 
V, (6) 0 


C<2 


C+ co 
Tellus XII (1960), 4 


INSTABILITY OF EKMAN FLOW AT LARGE TAYLOR NUMBER 


with a suitably chosen smooth transition region 
at ¢= whose width is of order (mx)° but 
much less than one. 

The method of solution is as follows: Since 
the coefficients of (24) are analytic there will 
be five linear independent and analytic solu- 
tions. Asymptotic representations in |mx|> 1 
for each of these may be found in the region 
o < ¢ <a and this will be referred to as the 
main boundary. These expansions will be 
seen to be singular at the end points, and one 
then looks for different representations valid 
in the vicinity of &=o (called the lower 
boundary solution), and a set valid for &> x 
(called the upper boundary region). These 
must then be joined or otherwise identified to 
satisfy the end point conditions. 


(a) Main boundary solutions 


Equation (24) has one particular solution ¢, 
and five complementary solutions {4,, $9, da, 


Pa $5}. For Vo(¢) #0 we may obtain two of 


these as asymptotic expansions of the form 


Substituting this in (24) and equating the 
coefficient of the highest power of (mx) to 
Zero gives 


Vo [Vogi (6) + Vogi (6)] =iV, 


fi 
&1 (6) = constant x Vy (6) —iV, (25) 
0 
or 
iV, [de 
> a iR V2 


d,— constant x V,(C) 
As |mxz|—ov, foro<ê<zx (25a) 
When £ is small compared to unity (but still 
in the main boundary) V,=e-t sin = ¢ and 


(25 a) becomes, d,> 2 db, > ¢. Note how- 


ever that ¢/ (¢) is singular at =o. 


Asymptotic forms for {2, Ps, Pa, ps} in 
the main boundary will be determined using 
the W.K.B. approach, in which it is asserted 
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that there are rapidly (in ¢) varying com- 
plementary solutions of the type: 


= exp [(mx)' S(E, mx)] (26) 


When the inhomogeneous term from (24) is 
removed and (26) substituted therein one gets 
the following non-linear equation for the 
unknown function S(¢, mx) : 


(iS) 217 (5) - V35)]= 
- (mx) a [10 8”(S’)® -6iV,S"S" - VS + 
+ ViVi] - (mx) [108 (5) + 155'(S”)2 - 
Ser a oad _ 
- (mx)-2 (SY + 48’) 


(27) 


S is now expanded in a power series in terms 
of the small parameter (mx). Let: 


S(E, mx) = Sy(C) + 0 + (mx)—hS, (0) + 


+ (mx)... (27a) 


Then the leading term S)(¢) may be obtained 
by setting the right hand side of (27) equal to 
zero, thereby giving: 


sd IS)? - iVo]? =0 (28) 


Eq. (28) is a fifth order algebraic equation in 
er whose solutions are: 
We ; 


So, = constant 


(a 
{So2 Sos, Soa» Sos} = + al As (¢) do + 


+ constant (29) 
The root Sj,=0 corresponds to the slowly 
varying solution ,, as may be seen more 
clearly by obtaining the next approximation 
Sn. The two roots +0'V 9/2 are degenerate 
roots of the quartic, which are resolved into 
four distinct roots in the next iteration of 
(27). When either of the two roots in (29) 
is substituted into the right hand side of (27) 
the term containing (mx) is found to 
vanish. If this were not the case, the first cor- 
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rection in (27a) would be O(mx) "4 and not 
O(mx)-"». Although formal substitution of 
(27a) into (27) then gives an intractable non- 
linear differential equation for the first approxi- 
mation $,/(£), we know that the four rapidly 
varying eigenfunctions have the form: 


{Po Pa Pas Ps} = 
= exp + [v/s (mx) f Vohde) exp + [S,(¢) + 
+ O (mx) "] (30) 


If we content ourselves with the zero order 
solution, then 


G 
& ~ constant x exp [i" (mx): [ V5 4€] 
0 
(31) 


The symbol (~) denotes equality in the sense 
that the fractional departure of the logarithm 
of & from the log of (31) approaches zero as 
| mz | > 00. 


i =2—k(1 +1) for mx >o 


(b) Upper boundary solutions (Vy = 0) 


We have assumed that V, is zero for € some- 
what greater than x. In this region the eigen- 
functions must satisfy: 


BY +4 =0 (32) 


which is obtained by setting V,=o in (24). Eq. 
(32) has two solutions which decrease ex- 
ponentially with height, one which is constant, 
and two which increase exponentially with ¢. 
The two of the main boundary solutions 
(presumably one from each of the degenerate 
solutions in (29)) which join to the exponenti- 
ally increasing upper boundary solutions must 
be rejected to satisfy the condition #’(®)=0; 
and the coefficients of the remaining functions 
must be chosen to satisfy: 


p (0) =2 
after the proper join is made at € =o. 


(c) Sub-boundary layer 


When the zero order roots (29) are sub- 
stituted into the right hand side of (27) and 
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~€+0 we note that the result becomes infinite. 


The dominant term is the one containing the 
fifth derivative of S. Using Vo(¢)=¢, one 
finds that as € approaches zero: 


(mx)-2 SY = - 2 (i) (mx)'2¢-h 


By comparing this with the left hand side of 
(27) as Ê +0, it can be shown that (29) is a 
valid asymptotic approximation in the range 
63> e3=|mxz|—1. In order to find a repre- 
sentation of the eigenfunctions for values of & 
which are not large compared with e we 
return to (24), stretch the vertical coordinate 
by the transformation €=7|mz|—*, and 
perform an asymptotic expansion by again 
letting |mx|>oo. With V,=|mx| "7 this 
transformation of (24) leads to the following 
fifth order equation: 


Er glad d 
= i up as N ve 


N 


3 
we introduce mxd = 

ini p=y 

and get the following equation for the sub- 


boundary layer: 


E et Peo d + 
dn? N dn? me Tin Jr- 


3 
Since ie = (mx) ce 


(34) 


Note that y,= -i and y, = constant x 7 are 
solutions to (34) which join on to ¢, and ¢, 
in the main boundary. More important, we 
must require that the solution of (34) joins (31). 
Using a W. K. B. technique similar to a) above 
it is readily shown that 


20 
y ~ constant x SIE ley | 
3 


vh=2h(r +i) 


(35) 
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is a logarithmically convergent asymptotic 
solution of (34) for n>1. On the other hand 
for e<€<1 eq. (31) becomes 


y ~ constant x exp E ple (mE | (36) 


which is identical to the sub-boundary solution 
(35). We can also obtain power series (in 7) 
representations for y which are valid for 7 <1. 


d) Piecing together 


In sections a), b), c) above we have obtained 
information about the form of the eigen- 
function except in the vicinity of the two 
transitions points 7 = 1, €=z. In order to join 
the functions in the various regions it would be 
necessary to obtain complete solutions of (34) 
over the whole range of 7 and to have the 
next approximation to the main boundary 
function. This we have been unable to do? 
and must therefore resort to an artifice which 
is consistent with the logarithmic accuracy of 
the above eigenfunctions. 

We ask the following question. Is it possible 
to find a single functional representation which 
is analytic at all values of ¢ and which: 

(i) satisfies the boundary conditions 

(ii) satisfies the differential equation in the 
main boundary by reducing to some linear 
combination of the zero order asymptotic 
solution in (31) fore<€ <a 

(iii) satisfies the upper boundary layer equa- 
tion (32) for &>n 

(iv) satisfies the lower boundary equation 
(33) for 7< 1 and forn> 1? 

Although such a fit may be expected to 
describe the velocities in the transition zones 
£ =eand £ = x with even less accuracy than in 
the main boundary, the size of these regions 
and the velocities in them are relatively small 


oo 
so that their contribution to /vd£ will be 
0 


negligible when |mx|> 1. 
If the two zero order solutions in (30) are 
combined one obtains the function cosh 


is 
[(mx) hi [ Vord£] which is much greater than 
0 


2 See RABENSTEIN (1958) for recent mathematical in- 
vestigations of these difficult joining and identification 
problems. 
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unity in the main boundary. Accordingly the 
function 


$ = Afı - cosh[(mx)#i# J Vo Hd (37) 


reduces to (31) and is an asymptotic solution 
in this region. In contrast with the exponential 
function (31), eq. (37) is analytic at =o 
because in the vicinity of this point 


ë 5 3 
TRA eels aoe 
LV dé 3° (: Bre...) 
and 


& 
cosh [ (mx): [ Ve dE] = 
0 


i /2\2 2 
=I +i(2) Ze nn +. 
PAS 10 


Therefore, if we choose A= Siecle then 
2imx 


(0) =2 and the function 


b= [1 - cosh {(mz) if Vehd) 
(38) 


satisfies the conditions (i), (ii) above. Further- 
more for > x eq. (38) approaches a constant 
and it therefore satisfies condition (iii). For 
e<€ <I eq. (38) reduces to (36) or (35), while 


for o<C<e it reduces to — Ci + 0(7*) 


which is one of the polynomial solutions of 
(34). The function in (38) therefore satisfies 
conditions (iv) above. 

An alternative interpretation of the sense in 
which (38) is an approximate solution is as 
follows: If it is substituted in the original 
differential equation (24) and the r.m.s. value of 
the resulting expression is divided by the r.m.s. 
value of one of the large terms in (24) (say 
dV), the result approaches zero as (mx#)-1. Thus 
the error is similar to the distortion which 
would be produced in the eigenfunction 
(satisfying the given boundary conditions) if 
an additional force were imposed on the 
physical system whose magnitude was small 
compared to the viscous force by a factor 
(mx) =!. 
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To get the approximate value of {vp dE we 
0 


set £ = co in (38) and obtain: 


31 Mae er 
~ —_[ı - cosh {(mx) i'l f V dE 
Boom) = [x - cosh (mn Ve nde] 
(39) 
This expression “converges” on the exact 
solution (for positive definite 1) only in the 
sense that the percentage deviation in the 


logarithm approaches zero as |mx|— oo, Or 
stated otherwise we have: 


dome 
31 1}. ot = 1 
nn [1 - cosh {(mx) hi ” Voledl + O (x)}] 
(39a) 


When the real and imaginary parts of (39) 
are separated one obtains: 


We note that (40) only depends on the mean- 
root value of the radial velocity profile and 
therefore is not very sensitive to detailed 
structure. The value of y which is used in the 
subsequent calculation was approximated by 


taking: 


TT 27 
y = fete (sind) dd + [ e=2l2 sin él dé 
0 n 
Vie Tide OR OA 


(42) 


where the last term is an error estimate. 
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6. Marginal Stability and the Associated 
Transports 


Eq. (39a) is an oscillatory function of |mx| 
and therefore Red (co, mx) is certainly positive 
for some large values of the parameter. In 
view of this, we conclude from the growth 
rate formula (18) that cyclonic and anticyclonic 
Ekman flows are unstable to body-boundary 
modes at sufficiently large values of the Taylor 
number and the Reynolds number. 

The asymptotic criterion for marginal sta- 
bility is given by (22), and using the ap- 


proximate relation (40) one then gets: 


I mx)? 
42 = —— li rm. (43) 
ver 
; 75 en 
for |mx|> ı 
or 
dre (xp)$ 
te (sinh zu) sin zu (44) 
A 
where zu=y = 


Fig. 2 is a plot of (44) showing the first few 
of the infinite set of discrete branches which 
are associated with this equation. The ordinate 
is proportional to the value of x=RE" 
which is necessary for the marginal stability of 
the mode with non-dimensional radial wave 
number m, and the abscissa is proportional to 
|mx|'». Shown below this is the numerical 
value of |mx| when y=1, which is approxi- 
mately the value of (42). The dashed part of 
the first branch at small mx indicates the failure 
of the asymptotic expansion, since the flow 
must be stable at mm=o. A third axis (E~/*) 
has been drawn to emphasize the fact that the 
condition of marginal stability is determined 
by a surface in x, m, E space (in addition to 
the degenerate wave number /) and that the 
curves which have been drawn represent the 
intersection of that surface with the plane 
E=oo, The heavy vertical lines xm = con- 
stant (in Fig. 2) are asymptotes to the various 
branches I, Il, Il],... and each unstable 
branch is separated from the next by stable 
regions in which the energy transfer is from 
the disturbance to the mean flow. If these 
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Figure 2. The section of the marginal stability surface 
for very large Taylor number showing parts of the four 
largest scales of motion. x = RE", y = 1.1, m= non- 
dimensional wave number. Only the minimum value for 
scale IV is indicated. The dashed part of I indicates its 
unreliability at values of |mzx| which are not large com- 
pared to unity. The inset labeled M.T. (amplitude 
arbitrary) is the vertical momentum transport in the 
interior of the fluid as a function of mx. The maximum 
point of this curve corresponds to the point which has 

been used in estimating an upper bound for R.. 


curves are plotted with m (rather than mx) 
as abscissa and x as ordinate then the asymptotes 
to the various branches will be the hyperbolae 


I : 2 
x=— x constant. Consider the stability, for 
m 


a given value of m, as x (i.e. U) is increased. 
A line of constant m will intersect each branch 
at two values of x, and the disturbance will 
grow when the value of x lies between these 
two. As x is increased further the energy 
transfer is reversed until the next branch is 
reached where there will be two more inter- 
sections, each corresponding to marginal 
stability. However the intersection of a constant 
m line with the (n +1)" branch occurs at larger 
values of x (and mx) than the corresponding 
intersection with the nt? branch, and therefore 
the associated eigenfunction will be a more 
rapidly oscillating function of ¢. This state 
of affairs is similar to the infinite set of discrete 
eigenfunctions which is found in the theory of 
thermal (Rayleigh) convection at marginal 
stability, and by analogy each branch in Fig. 2 
will be referred to as a ‘‘scale” of motion. Each 
scale has a preferred disturbance, which in the 
asymptotic limit E> oo is represented by the 
(relative) minimum point (x,,m,). In this 
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limit the minimum Reynolds number neces- 
sary for the onset of the nt" scale is given by 
R,= %,E" and the preferred wavelength 
in Centimeters is L,=27#,1E "4 H=27m;1 
E‘hhg. In contrast with the case of thermal 
convection, where the least rapidly (in z) 
varying eigenfunction always occurs before 
(at lower Rayleigh number) the higher order 
eigenfunctions, Fig. 2 indicates that the prefer- 
red scale is a function of the Taylor number. 


Since lim 4,41 < %, (with the possible ex- 
E>o 


ception of n=1) the (n+1)* scale should set 
in at a lower Reynolds number than the mp 
scale, at sufficiently large Taylor number. It 
is therefore necessary to obtain the functional 
dependence of x, on E when 1<E <(m!) 
(cf. remarks at end of sec. IV) in order to 
obtain a complete picture of the marginal 
stability. 

The inadequacies of Fig. 2 in the region of 
practical interest may also be seen from the 
following calculations. First consider scale II 
for which the minimum is x, = 13.6 at m,=5.8. 
According to our error estimate in the Taylor 
number expansion this value of x, can only be 
relied upon when E’/?~m,°~4 x 104, whereas 
a typical laboratory value is E'?~ 5 x 101. In 
scale I however the values of m are much 
closer to unity. For example the point u = .8 has 
an ordinate x = 11 and corresponds to mx =12.6 
when y=1.1. In this case m=1.1 and m® is 
only 1.7 which is negligibly small compared 
to E’'"=so. Although we cannot resolve a 
minimum point for scale I because of the failure 
of the second asymptotic expansion at small 
mx, it is believed that the value mx = 12.6 is 
sufficiently large to justify using that point as 
an estimated upper bound for the minimum x 
in scale I. Thus 


a = HIS 11 (45) 

Our discussion has been confined to the 
linear problem but it is of interest to compute 
the vertical transport of (angular) momentum 
by the Reynolds stress in the interior region 
of the fluid, in order to gain a picture of the 
tendency of the mean field modifications 
arising from the instability. Referring to Fig. 1 
it is seen that the downward transport of 


cyclonic x’) momentum is — w’u', where the 
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average is taken over one wavelength in both 
horizontal directions. For vertical distances 
much greater than the Ekman boundary 
thickness, the velocity components are found 
from (3) and (9) with u., = 1, and are given by: 


u = Ree® 


w = Re -iE"k[m [vyod& + 2 (mveo + IIe? 
0 


6 = At +iymE™! +ixlEh, 
Since w’(1, x, y) =0 we find by virtue of (16) 
w’ (I —2Z, x, y) = ReizE~smdg (oo, mx)e° 


where z2+Z2=1. 


Therefore the downward transport of cyclonic 
momentum at a distance Z < 1 beneath the free 
surface is 


x 


> El Me 
2 E -cosh( >) x 


Ya 
y )] e2t Rea 


in which the asymptotic expression (41) has 


mx [le 
been used. Therefore if cos y | — 
2 


(46) 


= COS nu 


is negative the downward transport of mo- 
mentum has the same sign as x. Referring to 
fig. 2 or (44) it will be seen that the minimum 
points for scales II, III... correspond to 
angles zu which are in the second quadrant. 
When x>0 (cyclonic rotation of the basic 
flow) the preferred infinitesimal disturbance 
transports cyclonic angular momentum down- 
wards from the interior regions and tends to 
produce a uniform reduction of the mean 
momentum here. Similarly, when x <o the 
Reynolds stress produces a downward trans- 
port of anti-cyclonic momentum. In either 
case the preferred disturbance tends to transfer 
momentum down the mean gradient. It should 
be pointed out however that in each scale 
there are less favored disturbances for which 


MELVIN E. STERN 


cos zu is positive aid. which transport mo- 
mentum up the mean gradient. In inset labeled 
M.T. in Fig. 2 is a plot of eq. (46) as a function 
of mx for the largest scale of motion. Although 
the amplitude is arbitrary its sign is positive 
and there is a maximum in the vicinity of the 
point which has been chosen to estimate an 
upper bound for x,. It is also of interest to 
compute the radial transport of momentum 
due to the correlation of w and v’ in the in- 
terior. The amplitude v is obtained from 


(16) and 1,9 = 1. One then computes u’v’ = 
a, Red( ) —21/m] ES pe 2) 
= 2Red(co, mx = 
f 1 Von 


at marginal stability. This is positive for axial 
symmetric disturbances or for 2/ < m5, and in 
these cases the growing disturbance tends to 
transport cyclonic angular momentum towards 
the center of rotation. 

We should also like to know to what extent 
each of the two components (Up, Vo) is 
responsible for driving the perturbation. This 
question is most readily answered from eq. 
(21), which is the boundary layer equation 
before the functional relation between U, and 
V, is introduced. Note that U, only appears 
in one term which subsequently vanishes by 
virtue of (23). However this term is also 
relatively small at large values of mx and 
would be neglected in sec. 5 on that score. 
It is therefore the momentum terms which are 
associated with V, which are of importance in 
the boundary layer at large values of mx, and 
consequently it is the work done by the Rey- 
nolds stress on this radial component of mean 
velocity which balances the dissipation in the 
boundary layer and in the interior. 


7. Summary of assumptions and conclusions 


The model abstractions and mathematical 
approximations which were used in this theory 
are listed below. 

i) The basic Ekman flow has no horizontal 
divergence and is only a function of the 
vertical coordinate. Such flows can only be 
approximated in the laboratory due to the 
circular geometry of the rotating vessel and 
the side walls. The influence of a basic gravita- 
tional stability has not been considered, al- 
though it is recognized that in the geophysical 
context this quantity may have a controlling 
influence on the scale height. 
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ii) The slope of the free surface which must 
accompany the rotation has been neglected in 
applying the boundary condition on the 
vertical velocity. This is a valid approximation 
if the rotational Froude number does not 
exceed unity. 

iii) The pivotal approximation in this study 
is the boundary layer separation of the eigen- 
mode and the subsequent expansion in the 
Taylor number. Since the associated error is 
of order m®E-"": this is quite adequate for the 
largest scale of motion (m~ 1), when E's Sr. 
However for the smaller scales (II, HI...) the 
preferred value of m is much greater than 
unity and E’: must therefore be much larger 
than m® in order that the asymptotic result 
be valid. 

iv) The boundary layer equation for the 
eigenfunction which is obtained from this ex- 
pansion (iii) is a function of only one param- 
eter (mx). Slowly “converging” approxima- 
tions for the necessary vertical integral of this 
function have been obtained using a W.K.B. 
approximation for the limiting case |mx|> 1. 
This is inadequate for computing growth rates 
at the long wave-length side of scale I, and 
does not give a resolution of the minimum 
critical x for the largest scale of motion. 

v) In connection with iv) the Ekman profile 
V,=e” sin € has been smoothed and the 
eigenvalue found for arbitrary positive definite 
V5. Although it is possible that additional 
modes may be associated with the neglected 
ripple component in V, it is believed that the 
minimum critical Reynolds number will only 
be affected to the extent of the r.m.s. departure 
of the realized profile and that positive definite 
function which was chosen to fit it. 

vi) We have only considered the growth of 
infinitesimal perturbations and leave unan- 
swered the following question. Is it possible 
that modes of finite amplitude produce a 
momentum transport which alters the mean 
field in such a fashion that the disturbance 
becomes unstable at lower Reynolds number 
than a corresponding one of infinitesimal 
amplitude? It is worth pointing out that the 
linear theory of Poiseuille instability, which is 
much more accurate than this calculation, does 
not give quantitative agreement with observed 
minimum critical Reynolds numbers. There is 
reason to believe that finite amplitude effects 
are responsible for this discrepancy. 
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vii) We do not know the stability charac- 
teristics of the pure boundary modes (i.e. 
modes of the first species which were referred 
to in the Introduction). If these occur first, 
then we would have to know the subsequent 
modification of y to determine the onset con- 
dition for the larger scale body-boundary 
modes (assuming a separability of the two 
stability problems on the basis of the difference 
in scales). 

The following conclusions have been drawn 
from the theory: 

i) A non-convergent Ekman boundary flow 
is unstable to body-boundary modes at suffi- 
ciently large values of the Taylor and Reynolds 
number. These draw their energy primarily 
from the cross-isobar component of the mean 
field and may be expected to produce a direct 
modification of the interior geostrophic flow 
depending on the difference of the Reynolds 
number and the critical R. 

ii) The estimated upper bound for the mini- 
mum. critical Reynolds number is given by 


Rie lene TT 


Thus for a typical laboratory value of E = 2,500 
the Ekman flow should not be stable at 
Reynolds numbers as large as 80. It is con- 
ceivable that instability may set in at much 
lower values for reasons mentioned previously. 
ii) The radial wavelength of the preferred 
disturbance which is associated with ii) is 


Ly= 20m, E“hs 


where m, is a numerical constant which can 
only be estimated as m,=10°. The azimuthal 
dimension of this disturbance should be much 
larger than the scale depth. 

iv) At , the disturbance with wave number 
m, is marginally stable while disturbances with 
slightly different wave numbers are damped. 
There are also (relative) minimum Reynolds 
numbers corresponding to the marginal sta- 
bility of the higher order eigenfunctions. The 
preferred cigenfunction, or the smallest of the 
set of relative minimum Reynolds number, 
depends on the Taylor number. At sufficiently 
large E the (n +1)" scale is preferred over the 
nth scale. There may be certain critical values 
of E where two such disturbances, with 


414 


markedly different wavelengths, can be ob- 
served at the same minimum critical Reynolds 
number. 

v) The relative magnitudes of the perturba- 
tion velocities in different regions may vary 
considerably, depending on the values of the 
Taylor and Reynolds numbers. The radial 
component near the free surface is smaller 
(by a factor of E~) than the velocities in the 
boundary layer, while the latter are large com- 
pared to the interior azimuthal velocities by 


order |#(®, xm) |. 


8. Experimental Evidence 


The author’s interest in the preceding sta- 
bility problem came from an earlier experi- 
mental study of the transient Ekman flow 
which is produced when a shallow cylindrical 
vessel of water is suddenly decelerated from 
one rotation rate to another. If the initial 
equilibrium rate w is decreased by a small 
amount Aa, the fluid will rotate at a uniformly 
greater rate than the container, except at the 
bottom (and the rim) where a quasi-steady 
and convergent Ekman layer will be formed. 
Above this, for reasons of mass continuity, 
there is a weaker divergent radial flow which 
is responsible for the local decrease of angular 
momentum in the interior. If rotationally 
stable floats be placed on the free surface, their 
total relative angular deflection (©), during the 
readjustment process, is an inverse measure of 
the mean decay rate. These experimental 
quantities will be compared with values pre- 
dicted from the Navier-Stokes equations on 
the basis of the assumption that the flow is 
always axi-symmetric and laminar. We would 
like to ascribe systematic departures in this 
comparison to an asymmetric instability, even 
though it is recognized that there are con- 
ceptual difficulties in defining stability for 
inherently transient mean fields. 

The cylindrical tank which was used in these 
measurements had a radius (A) of 52 cm, and 
was centered and leveled on a rotating platform 
built by von Arx (1952) for ‘oceanographic 
model experiments. It was driven by a variable 
speed motor in the range Aw/a= 1), to 1}, 
with most measurements in the vicinity of 
Aw /w@ =1/,, and @=1 rad/sec. Fluctuations in 
motor speed of 1% were noted and this 
produces a 10 % uncertainty in the parameter 
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Aw/w at average conditions. However the 
period of these fluctuations is much smaller 
than the four—five minutes which occupy the 
entire decay process and therefore these fluctua- 
tions will be partially averaged out if one uses 
four—five minute averages of the initial and 
final motor speeds to compute Aw. This is 
believed to be the limiting factor in the un- 
certainty of the experimental data which is 
estimated to be less than 5 %. The bottom of 
the container was fitted with a glass plate and 
covered to a nominal depth of H=5 cm 
(kept constant in all runs) with tap water. The 
temperature of the latter was only measured 
at the start and finish of a day’s experiment and 
the viscosity was computed from the mean. 
Since the dissipation of total angular mo- 
mentum depends on the square root of the 
viscosity the uncertainty in our measurements 
due to this factor is o—3 %, less than or equal 
that due to the fluctuation in motor speed. 
After some trial, it was decided to use four 
1/,’’ circular cardboard floats as an indicator of 
the relative angular deflection of the water 
near the free surface. These were tested before 
the experiment with regard to orbital sta- 
bility and relative deflection, when the fluid 
and container were in solid rotation. In a 
time interval (4—s5 minutes) corresponding 
to the length of each experimental measure- 
ment of ©, no significant angular deflection 
could be detected. The four floats were initially 
set 90° apart, and their distance from the center 
of rotation ranged from A/4 to 34/4 with 
most measurements near mid-radius. A glass 
cover plate placed on top of the container 
isolated the surface of the fluid from the drag 
of the air in the room, and a grid ruled upon it 
served to measure the initial and fina] position 
of the floats. The average relative angular dis- 
placement was one revolution and the un- 
certainty in obtaining the coordinates is negli- 
gible. 

The data so obtained was plotted in Figs. 3 
and 4. The ordinate in these figures is es- 
sentially the total angular deflection per unit 
change of angular velocity, since for small 
values of Aw/w the logarithmic ordinate be- 


Aw 


—1 
comes -0() . In this limit the abscissa 


2 A2 
reduces to E'» (1+F/4), where F = (2 = ) is 
8 


Tellus XII (1960), 4 


INSTABILITY OF EKMAN FLOW AT LARGE TAYLOR NUMBER 


the rotational Froude number. It therefore is 
essentially a measure of the mean rotation rate 
since H was held constant in the experiments. 
It is easy to show by integrating the linearized 
Navier-Stokes equations with respect to time 
that the relationship between these two limiting 
quantities is linear to the extent that Aw/w, 
Eh} and H/A are much less than unity. Since 
Aw could not be controlled by the nominal 
motor setting but had to be measured after the 
fact, we have extended this calculation to 
include larger values of Aw/w. The elaborate 
coordinates in Fig. 3 and Fig. 4 are an attempt 
to reduce the data so that they fall on a single 
curve, for all Aw/o. If the relative motion were 
laminar, then our points should all fall on the 
heavy straight line which, under present 
operating conditions, differs by no more than 
1 % from an exact initial value solution of the 
Navier-Stokes equation. This calculation has 
not been included herein because it is long and 
not especially pertinent. It is based on a bound- 
ary layer expansion in E~‘ and includes finite 
amplitude effects, correct to order (Aw/w)ÿ, as 
well as the variation in the height of the free 
surface and the vertical velocities associated 
with the convergent flow. 

Each plotted point in Figs. 3 and 4 is the 
average displacement of the four floats in a 
given run and the thin vertical line gives the 
extremes. The geometrical symbolism used in 
plotting the points reflects variations in R 


(muipty OO pr by A/H= 10 to get R) 


For relatively small values of E and all values of 
R the experimental points are, within observa- 
tional error, in agreement with the straight line 
predicted by the symmetric theory. However, 
Fig. 3 shows a systematic departure from this 
line beginning at values of E' of about 50 
(the maximum value of the rotational Froude 
number in these experiments was 1.6). All 
values of the Reynolds number are associated 
with a larger decay rate than one would expect 
if the motion were laminar, and there is less 
reproducibility and coherency of the points 
corresponding to fixed R. This increase in the 
average rate of dissipation is attributed to an 
“instability” in the “slow” decay of the laminar 
Ekman profile. The peculiar aspect of the 
transition zone (abscissa s0—s5) in Fig. 3 is 
that the critical parameter appears to be the 
Tellus XII (1960), 4 


Figure 3. The average rate of decay of an Ekman flow 
produced by decreasing the rotation rate (w) of a shallow 
cylinder by the amount Am-@ is the total relative 


angular deflection (in revolutions) of floats placed on 
2 A2 
the free surface. F= ws 


is the rotational Froude num- 

ber. To get maximum boundary layer Reynolds number 

multiply figures in legend [Sex| by A/H= to. Heavy 
o 


straight line is theoretical computation of the decay 
assuming laminar mean motion. Thin vertical line 
through the experimental points indicates the extremes 
in the displacement of four surface floats and the broken 
curve is an envelope of the data (In = log e). 


Taylor number rather than the Reynolds 
number. The latter range from values of 
40—130 in a small interval about E’:= so. 
In Fig. 4 (anti-cyclonic relative motion) the 
data are in satisfactory agreement with the 
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Figure 4. Same as Fig. 3 except Aw>o. The relative 
motion is now anti-cyclonic and the Ekman boundary 
flow diverges from the center of rotation. 
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laminar theory to somewhat larger Taylor 
numbers. Above abscissa of 60 the dashed 
envelope of the experimental points indicates 
a departure from laminar decay but the data are 
inadequate to infer any systematic behavior. 

In addition to these measurements, sepa- 
rate visual observations were also made to 
detect the instability. A thin layer of per- 
manganate dye was diffused on the bottom of 
the rotating tank and concentrated as an annular 
ring near the circumference. When the rotation 
rate is suddenly decreased, this ring of dye con- 
verges towards the center of rotation with the 
Ekman boundary flow. During the first two 
revolutions after the change in motor speed 
the inside of the converging permanganate 
front retains its initial circular form. At this 
time, however, the circular front (located at 


radius ‘e to A/2) begins to fold along about 


one-fourth of the circumference, as if there 
were an additional surge of fluid towards the 
center. Seconds later, roll striations appear with 
axis parallel to the fold, as was indicated by 
areas of varying dye concentration. These rolls 
then spread from this area and one may see 
the striations in many regions of the fluid in 
complex (but rather stationary) orientations of 
the axes. Considerable variation in the size of the 
small dimension of the roll have been detected 
in photographs (not included herein). These 
range in size from less than to much less than 
a centimeter while the large dimension of the 
roll has a wavelength of perhaps one-third 
the circumference. It has been reproduced 
many times? even in the region of Fig. 3 below 
abscissa of 50! I do not believe that the dye 
concentration was sufficient to introduce any 
extraneous effects in these qualitative observa- 
tions. No systematic asymmetry was detected 
in the surface motion! 

For the experimental conditions correspond- 
ing to Fig. 3, the following conclusions have 
been drawn. The quasi-steady Ekman motion 
which is produced by decelerations in the ro- 
tation rate of the container is “unstable” for 
certain values of the parameters. The dye ob- 
servations indicate that the energy source for 


# Similar disturbances have been seen in the large 
(20 foot) tank at Woods Hole. The depth of water was 
greater than in the above case while the basic rotation 
was about 1/, rad/sec. 
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the disturbance is in the boundary layer of the 
mean flow. However the observed transition 
point (abscissa of so in Fig. 3) is much more 
sensitive to changes in w than Aw (or U), and 
therefore it cannot be entirely represented in 
terms of a boundary layer Reynolds number 
R=U(rw) "h. It is suggested that the Taylor 
number is an equally important non-dimen- 
sional parameter, and that the scale depth plays 
a role in the instability. The independent visual 
observations also suggest that the flow was 
unstable at abscissa less than so. Therefore the 
transition which is detected by the quantitative 
measurement may mark the onset of a new 
scale of disturbance, which has a relatively 
larger influence on the mean field at these 
operating conditions. The dashed line, which 
is the envelope of the experimental data in 
Fig. 3, suggests a tendency towards stabiliza- 
tion as the value of the abscissa increases to 
too. It would be interesting to see if a new 
transition is observed as the Taylor number 
is increased. Regarding the data for the anti- 
cyclonic flow (Fig. 4), I do not feel it is ade- 
quate to attempt to draw conclusions at large 
E’: concerning the proposed instability. 

It should be emphasized again that there are 
substantial differences between the theoretical 
model and the laboratory experiment. For 
reasons which have been mentioned a quantita- 
tive comparison of results cannot be made at 
this time. Whether or not there is any intimate 
connection in mechanism, I think the two 
studies are complementary in suggesting rather 
unexpected couplings between the friction 
boundary and the quasi-geostrophic flow in 
the interior of a rotating fluid. 

Prof. A. ARONS and his students at Amherst 
College (1961, in press) have performed an 
experiment in which an Ekman flow is formed 
by forcing water into a rotating tank at the 
inner rim and letting the water level rise. They 
have informed me of axial symmetric instabili- 
ties which have a structural resemblance to 
those which are discussed in this theory and 
whose radial wavelength is proportional to 
hgE'h. In the range of E from 1.5 x 10% to 
2 x 10* they find a constant of proportionality 
m=3.0+20%. However, these instabilities 
are observed at critical R ranging from 1.5 
to 3.4, considerably below the upper bound 
which has been estimated by this theory. 

Dr. A. Faller of the Woods Hole Oceano- 
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graphic Institution is also experimenting with 
Ekman instability, in a geometry more like the 
one to which this theory pertains. I wish to 
acknowledge with gratitude the discussions in 
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which his ideas and preliminary results were 
freely given. The work reported here has been 
supported by the Office of Naval Research 
(U.S.) under Contract Nonr 2196(00). 
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Abstract 


The equations of motion are developed in a form appropriate to flow parallel to a uniformly 
sloping bottom. Some cases of two-layer flow are considered in which the flow is concentrated 
in the lower layer, as in the overflow over a submarine ridge. Whereas in the absence of friction 
such currents would tend to follow the contours of the bottom, the effect of bottom friction 
is to deflect the current to the downslope side of the contours. Numerical solutions show that 
in some cases the angle of deflection might be as large as 30°, but in most circumstances it 1s 
probably less, of the order of 5° or 10°. An indication is given of how the theory may be 
applied to the case of density varying continuously with depth, as in a real ocean. Some general 
conclusions, which are relevant to the interpretation of observational data, are drawn. 


I. Introduction 


The flow of water over a submarine ridge 
into an oceanic basin is frequently the source 
of supply of an important water mass. The 
overflowing currents differ from other oceanic 
currents in that they appear to reach their 
maximum velocity at, or near to the bottom, 
the overlying water being almost at rest. 
The flow of such currents is likely to be in- 
fluenced to an unusual extent by the contours 
of the bottom and the effect of bottom friction. 
Observations indicating the presence of such 
flow into the North Atlantic over the Iceland- 
Faeroe and Greenland-Iceland Ridges have 
been described by DIETRICH (1956, 1957). 
STEELE (1959) has also described observations 
of deep overflow across the Iceland-Faeroe 
ridge. The process of cascading (Cooper 
and VAUX, 1949) would also seem to imply 
the existence of currents of this type. 

The equations of motion will be developed 
in a form appropriate to flow parallel to a uni- 
formly sloping bottom. Solutions will be 
given for some cases of two-layer flow, the 
water in each layer being of uniform density 


and the thickness of the lower layer small 
compared with the total depth of water. The 
net flow in the upper layer will be assumed 
to be zero. An indication will be given of 
how the method may be applied to the case 
of density varying continuously with depth, 
as in a real ocean. 


2. Equation of motions 


Consider, in the first instance, flow on the 
side of a ridge with a straight horizontal crest 
and a downslope at a constant inclination « to 
the horizontal. Let rectangular axes XYZ be 
taken, with OY along the crest of the ridge, 
OX down the slope, at an angle « to the hori- 
zontal, OZ perpendicular to OX and OY 
and directed upwards, as in Fig. 1. Let OY 
be orientated at an angle © to the west of 
north. 


Let u, v, w be the velocity components, 
p pressure, © density, ¢ latitude, © angular 
velocity of the earth’s rotation, X, Y and Z 
the components of external force, other than 
gravity, per unit mass. 
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Fig. I. System of co-ordinates: OY is parallel to crest of 
ridge; OX is parallel to line of greatest slope. 


Me ae nu? +v—+w 

m “Ot ds dy 9 
Then the equations of motion may be written 
in the form 


EM ee v (cos asing + 
— = N 
Dt 


+ sin «sin © cos ¢) + wcos © cos d | 


= ss + gsina+ X 
0 0x 
2 20 I (cos&sind + sin «sinO cos) + 
Dt | (1) 


+ w (sin «sind — cos «sin@ cos 9) = 


rn 
ey 
Fi = 20 Jucos0 cos b+ v (sin «sind — 
0 
- cosasin 8 cos )| = =. foe Zi 


with the equation of continuity 


Ou dv dw 


ae ah deal (2) 
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In considering flow near the bottom, w will 
be small and terms in w may be neglected. 


If, also, the inclination « is small enough (of 
order 10-?, say) for terms in sin « on the left 
hand side of (1) to be neglected compared 
with those in cos «, and if cos « is taken as r, 
equations (1) reduce to 


Du 1 Op ; ] 
Du | 
Er el 
Dt u ayn (3) 
— 20 cos p (u cos ® - vsin@) = See | 
-gcosa+Z | 


where f =2@ sin ¢. 


If the terms on the left of the third equation 
of (3) and the external force Z are negligible 
compared with the main gravity component g 
cos &, this equation reduces to 


Op 
32 ~ SQCOs a (4) 


corresponding to the normal hydrostatic 
equation. 

Allowing for shearing stresses tangential to 
planes perpendicular to OZ but neglecting the 
lateral shearing stresses, the force components 
may be written 


BAR y TOR, 


% @ de’ 


where F,, and F,, are the components of 
shearing stress by which water above a plane 
z=const. acts on the water below that plane. 


3. Equations for two-layer flow 


Consider an upper and a lower layer of 
density ©, and 6, respectively, the thickness 
of the lower layer being (x, y). Let the free 
surface be given by 


z=h)txtana+C 


as in Fig. 2. Let p,= atmospheric pressure. 
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Fig. 2. Co-ordinates for two-layer flow: 7 = thickness 
of bottom layer. 


Then, from (4), the pressure is given, in 


the upper layer (z >%), by 
Pi =Pat go, cosa(hy+xtana+f-z) (6) 


and in the lower layer (z < ) by 


p=pa+gcosa {e,(hy+xtane+l — 7) + 
02m - z)} (7) 
From (6) and (5) in (3), and integrating 


throughout the depth of the upper layer, 
it is found that 


Du, =e) atop; | 
D: 1= es ne 
Fox — Enz 
o,h 
ls 111 (8) 
REN rp RE Gg 3 
Di fu, an Eee at 
Fey = Fay 
CLR 
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where h,=h)+x tane-7 


n+hı 
Er 
Uy = | Hide ic. 
1 
n 


F,, and F,, are the components of surface 
stress and Fyx, Fny those of stress at the interface. 


It will be assumed now that p, is uniform, 
the surface stress zero and that the mean 
motion in the upper layer vanishes. Then (8) 
reduces to 


cos & 2 Enz 
4 ere 
= ax Of 
(9) 
gCOs a a = — Fry 
oy QUA 


In the lower layer, from (7) and (s) in (3) 
and then integrating throughout the depth of 
the layer and using condition (9), 


where g’= (- or u == [ude etc., and 
Q2 1 


Fox, Foy, are the components of stress at the 
bottom. 


If h,>%, the frictional terms in (10) reduce 
to 


(Ey = Fox)/Q2 and (Fry 27 Foy)[03N- 


: d 
Since w =0 at z=o and w = at z =, the 


equation of continuity for the lower layer 
gives 


ENT Re 
Fin PA ee N 
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4. Form of the frictional terms 


Let the resultant mean velocity in the lower 
layer be V, ice. 


V2 = u? + v2 


and let the combined shearing stress at the in- 
terface and at the bottom be related to V by a 
quadratic law, so that if F is the resultant fric- 
tional stress, 


F=ko,V? 


If it is assumed also that the stress F is in a 
direction directly opposed to the mean velocity 
V, then the stress components on the lower 
layer are given by 


Fox — Fnx = ko, Vu 
Foy = Fyy = ko, Vv 


Thus in (10), with h,>n, 


SIT I: kVu 
ais = (sina — cosa) En 


There is little experimental evidence to in- 
dicate the magnitude of k in these conditions, 
i.e. a bottom current flowing below water 
sensibly at rest. From a number of observations 
of tidal currents in comparatively shallow 
water, values of k in the range 1.5 to 4 x 107? 
and averaging 2.5 x10~* have been derived 
(summarised data have been given by Proup- 
MAN, 1953, and BOWDEN, 1956). From wind- 
induced currents in the Straits of Dover, 
superimposed on tidal currents, BOWDEN (1956) 
found k=3.8 x 10-3 and ROssITER (1958) 
found k=3.5 to 4.5 x107-%. WVYRTKI (1956) 
has summarised a number of observations in 
which the friction has been estimated from 
wind drift in the surface layer in the open 
ocean. The corresponding values of k are 
mostly in the range 1.0 to 2.7 x 107%, apart 
from some rather doubtful values of 10 to 
12 x10-3 for the Atlantic. 

In dealing with the flow through straits, DE- 
FANT (1955) proposed the value k=3 x 107? 
and found that this led to estimates of velocity 
of the correct order of magnitude for flow 
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through passages such as the Dardanelles, 
Bosporus and Straits of Gibraltar. The flow 
in these cases differs from the bottom currents 
envisaged above in that the current in the 
bottom layer is bounded above by a current 
flowing with comparable velocity in the oppo- 
site direction. The resistance to the bottom 
current is likely to be greater, therefore, than 
if the overlying water were at rest or if the 
current extended to the free surface. The 
irregular bottom topography in the straits 
probably increases the resistance further. 

DIETRICH (1956) has used the value k= 
3 x10-? in estimating the rate of flow down 
the southwestern slope of the Iceland-Faeroe 
Ridge, from density data obtained in the 
“Anton Dohrn” cruises. The flow on the 
Iceland-Faeroe Ridge has also been considered 
by Steele (1960), who has suggested that the 
high apparent value of k may be due mainly 
to the loss of energy from the bottom layer, 
associated with turbulent mixing of the bottom 
water with that in the layer above. 

At the other extreme, STONELEY (1957), 
in dealing with turbidity currents, has suggested 
a value as low as k=2 x10~ for a turbidity 
current with velocity of the order of 28 m/s. 
This value is based on an extrapolation of data 
obtained by Keulegan, in model-scale experi- 
ments on the flow of saline wedges, in which 
he found that k decreased with increasing 
Reynolds number. 


5. Solutions for steady motion 


If all accelerations relative to earth are 
negligibly small compared with the other 
terms in (12), one may put 
Du Dv 
Baumes: and then 

Sal on kVu 
— fr =g (sina - od) Er 
(13) 


a: ‚on kVv 
U= — COS & 
f £5 r 
on 


Putting = =0 in (11) 


SON On ou 2) 
ET fo) 
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(a) Layer of uniform thickness 
If the thickness of the bottom layer is uni- 


form, i.e. = 5 =0 (13) becomes 
-fv NA 
N 
(15) 
kVv 


Let tany= -u/v, so that 
= —Vcosy. 

Then, from (15), 
(16) 


F kV 
Area 
7 


and 


RV2\1 o’sina 
Ale oe 7 


Thus V may be calculated, given 9’, sin «, 
f, n and k and hence, also, tany, u and ». 
From (16) and (17), 


5 re: +o) 
Panel Mn 


me g sina ( oi) 
p= - 1% 
Up Jon 

Since # and v, as given by (18), are inde- 
pendent of x and y, equation (14) is satisfied. 
Thus flow in a layer of uniform thickness, 
with uniform velocity given by (18), is a 
possible form of motion. 

If conditions are uniform parallel to the 
crest of the ridge, ie. An/oy=o, dv/dy=o, 
and if Ju/dx is negligibly small, it follows 
from (14) that 0n/0x =o also, so that the prob- 
lem reduces to that just considered. 

The case of flow in a bottom layer of uni- 
form thickness has been treated by Steele 
(1960), who obtained a solution essentially 
the same as that given above. 

In the frictionless case, k=o, the above 
solution becomes 


(18) 


@ sing 
Ni 


The flow is then purely geostrophic and 
parallel to the crest of the ridge. 


Min 


= 0; 
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If, on the other hand, the friction is very 
great, so that 


HA then u>V, v>o 


fy 
Li ayo DT TAT 
u=V ( k ) 


and 


so that the flow is directed down the slope 
and its velocity is limited by the friction 
coefficient k. 


(b) No transverse component of flow 


If there is no flow parallel to the crest of 
the ridge, ic. v=o, u=V, and if du/dx is 
negligibly small, then from (14), An/ox=o. 

Hence in (13) 


ku? 


o=g sina -—— 
(19) 


= Mm 
re 2 


Thus 
ia Fe 
(ue 69) 
and # 
ce cé 
dy £ 


Thus, in order to maintain a flow directly 
down the slope, a transverse gradient, given 


by (21), is required. 


(c) More general case 


If the slope is small enough for the approxi- 
mations discussed in § 2 to be valid, the axes 
OX and OY may be taken in any two perpen- 
dicular directions in the plane of the bottom. 
In this case, let OX and OY be inclined down- 
wards at angles « and ß respectively to the 
horizontal. Then equation (13) would be 


replaced by 


yey ig RE 
— fv = gli, - — 
a N 
(22) 
an yi teh ee 
Che ee, 
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A 2) net 2 
where i,=sin« =. cosa, i,=sinß oe cos fp. 
x y 
i, and i, represent the components of slope of 
the upper surface of the bottom layer in the 


OX and OY directions. 


If the OX axis is taken in the direction of the 
resultant current in the bottom layer, then 1 = 
V, v = o and (22) becomes 


SRE 

eee (23A) 

[V=giy (23B) 
(14) becomes 

oe = 

5, (nV) =0 (23C) 


Thus V can be computed from g’ and the slope 
i, of the interface in the OY direction only. 
Then, if i, and 7 are known, (23A) enables an 
estimate of k to be made. 


From (23A) and (23B) 
tany == =— 24 
Y fn ( ) 


This equation, corresponding to (16) for the 
case of a layer of uniform thickness, gives the 
angle y between the current vector and the 
contours of the upper surface of the bottom 


layer. 


6. Numerical solutions 


To obtain an order of magnitude for the 
quantities involved, consider the values given 
by DierricH (1956) for section IIIb across the 
centre of the Iceland-Faeroe Ridge. In the 
above notation, for a position at a depth of 
700 m on the south-west slope, 0, =1.0277, 
09 =1.0280, sin «=1/125, 7 =50 m. Then g’= 
0.29 cm/s?, and taking ¢= 63°, f=1.30 x 1074 
os 

Considering solution 5(a), i.e. the bottom 
layer of uniform thickness, and withk =3 x 107, 
equations (16) and (17) give 

V=14.5 cm/s, tany =0.67, y =33.8° 
so that 
4 =8.1 cm/s, v=12.0 cm/s. 
If k=3 x107°, 
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then 


V=17.5 cm/s, tany=0.081, y =4.6°, 
U=1.4 cm/s, v=17.5 cm/s. 


In the purely geostrophic case, ice. k =o, 
V=17.6 cm/s, tany =o. 


Even with the higher value of k, 3 x 1072, the 
component of flow parallel to the ridge is 
greater than the component down the slope. 

Considering solution 5 (b), ic. with no 
transverse component of flow, and taking k= 
3 X107?, equation (20) gives u=19.6 cm/s 


7 
and then (21) gives = — 0.89 x 10-2. With 
4 
k=3 x10-3, (20) gives u=62 cm/s, and from 
(21), à 292 10. 


Thus the flow can be directed down the 
slope only if the thickness of the layer increases 
to the right of the flow with a gradient of 
order 1077, ic. 10 m per km. This could be 
maintained in a narrow channel or canyon 
but not on an open slope. A comparison of 
Dietrich’s sections IIIb and IIa does indicate 
an increase in the thickness of the bottom 
layer from southeast to northwest, but only 
of the order of 100 m in 100 km, ie. 1073, 
which is too small by a factor of 10 to maintain 
a purely downslope flow. It appears, therefore, 
that the interpretation of a section such as IIIb 
should be based on solution (a) rather than 
solution (b), with the major component of 
flow parallel to the contours of the ridge. 

The conditions postulated in solution (c), 
equations (23 A) and (23 B), appear to resemble 
those reported by DIETRICH (1957) in the 
bottom current flowing southwards at the 
foot of the continental slope off East Green- 
land. From the velocities computed by Dietrich 
across sections off Cape Dan and Cape Fare- 
well, the mean velocities in the bottom layer 
may be taken as 10 cm/s and § cm/s respec- 
tively. The increase in depth of the isosteric 
surface, 0,= 27.9, from the C. Dan to the C. 
Farewell section is about 500 m in soo km, 
i.e. 1073. Taking Ag,=0:1 gives 9° =0.095 
cm/s?, and putting i,=10-°, 4=7.5 cm/s in 
(23 A), leads to 


et 70 cmEl. 
n 
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n may be estimated at 200 m, corresponding 
to k=3.4x10 2. 

Thus if the flow in this current may be 
assumed to be in a steady state, it is consistent 
with a value of k of the order 3 x10~?. Then, 
in the C. Dan section, taking #—10 cm/s,n = 
250 m gives tan,y=0.105, y—6°. In the C. 
Farewell section, with u=s cm/s, 4=150 m 
the values tan y =0.087, y =5° are obtained. 

It is of some interest to compare the features 
of the bottom currents considered here with 
those of a turbidity current. Taking the 
values 


01=1.0, @,=1.4, corresponding to 9’ =280 
cm/s?, k=3 x10-°, n=100 m, V=28' m/s, 
J=1.1 % 10% 574, then tal y=7.65 y= 8255. 


If the current is constrained to flow down 
the line of greatest slope, the transverse gradient 
of its upper boundary, from (21), is dn/dy= 
1-1X10-%, which is about 1/10 of that as- 
sociated with a bottom current of the type 
considered above. 


7. Density varying continuously with depth 


Let it be assumed that there is some surface, 
at a distance 2’ from the bottom, on which 
the velocity and the shearing stress may be 
taken to be zero, as in Fig. 3. Then from equa- 
tion (3), at 2=2z, 


2 


Fig. 3. Co-ordinates for density varying with depth: 
velocity and shearing stress zero on the surface z = 2’. 
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Be 0 
de = sine >= 


(25) 


=> 


If p.’=pressure at the point (x, z’), then the 
pressure at (x, z) is given by 


p=pz +gcosap, (2’ - 2) 


where 0, = = ei 


ES 
2 


Hence, using (25), 


902 
de = gina + g cosa(z’ -2) 
Then in (3), 
Du I doz 
Dr -fr= —gcosa (z ie 
hele 
@ 92 
Ps 2 (26) 
= 2 er 
= + fu gcosa(z’-z) = 
Fr 10F., 
G da 


If the accelerations relative to earth are zero, 
integrating from the bottom to z=2’ gives, 


2! ‚ 


Dan ch £ COS & : 9027 a 
fi ss 4: 0 0x ree 


c 


0 0 
Fox 
of 
ae zt {(27) 
7 a0 
dam BENE fi 
i} x Dar 
0 0 


— 


of 


since For = Fy,= 0. 


The terms on the left hand side of (27) rep- 
resent the components of volume transport 
by the bottom current. The first terms on the 
right hand side represent the components of 
transport in the absence of friction, while the 
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second terms may be regarded as correction 
terms due to friction. Then (27) may be written 


BETT, 
Il. til, 


where Ty [ude, T,= | vde, 
0 0 


= 


(28) 


ee 
Ton ey] BOs (oh aida 


f Jey 
0 
ASIEN 
Toy = 7 [exe z)dz, 

0 
AT, =~ 

are 

(29) 

AT, = 7% 


of 

Thus the transport correction terms depend 
only on the components of bottom friction 
and the Coriolis parameter. The numerical 
examples discussed in $ 6 suggest that in many 
cases the correction terms will not exceed 
about 10 % of the transports computed neg- 
lecting friction. 

A procedure for the analysis of observational 
data may now be outlined: — 

(i) Assuming that the position of a surface 
of no motion parallel to the bottom may be 
estimated, the velocity compone.ts across 
two mutually perpendicular sections, each 
perpendicular to the bottom, may be com- 
puted, neglecting friction. This may be done, 
for example, by Helland-Hansen’s method 
and the terms Tox, To, determined. 

(ii) From the computed velocity at the 
bottom, the frictional terms Fox, Foy may be 
estimated, assuming a value of the coefficient k. 
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(ii) The transport correction terms AT,, 
AT, may then be computed from (29). 


8. Conclusions 


The following conclusions, relevant to the 
interpretation of observational data, may be 
drawn. 


(i) In bottom currents, flowing on a sub- 
marine slope, the effects of bottom friction are 
likely to be appreciable and should not be 
neglected. At the same time, the geostrophic 
effects must also be taken into account. 


(ii) In the absence of friction, the currents 
would tend to follow the contours of the 
bottom. The effect of friction is to deflect 
the current to the down-slope side of the 
contours. The angle between the current 
vector and the contours might be as great 
as 30° in some cases, but it seems probable 
that in most circumstances it will be less, of 
the order of 5° or 10°. 


(iii) It is unreliable, under these conditions, 
to attempt to estimate the velocity of flow 
from data in a single vertical section. Observa- 
tions in two sections at right angles should 
be obtained. 


(iv) There is at present considerable uncer- 
tainty in the value of the friction coefficient k. 
It is desirable, therefore, that in as many 
cases as possible, direct measurements of 
bottom currents should be made at the same 
time as the density data are being obtained. If 
reliable estimates of k could be established, it 
might then be possible to compute bottom 
currents satisfactorily from density data alone. 
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Abstract 


An iteration method has been developed which uses finite differences for calculating water 
level elevation and water transport. To this end, the equations of motion were integrated over 
the depth and transformed into an implicit difference equation system. The equations were 
solved by the aid of an electronic computer. Two main cases were calculated, and in one in- 
stance the water levels were checked with the real ones. 


1. Introduction 


When the hydrodynamic equations of mo- 
tion of liquids cannot be solved by analytical 
means, a numerical method can be applied 
for the calculation of specific results. Most 
scientists have avoided this method because of 
the rather burdensome numerical work in- 
volved. Nevertheless, such calculations can 
often be carried out with the aid of electronic 
computers. Numerical calculations relating to 
a sea-level problem in the North Sea were 
made by Hansen (1956) and FISCHER (1959). 
They took the vertically integrated equations 
of motion, and transferred these two-dimen- 
sion equations into difference equations. These 
equations were then used directly for calcula- 
tions of the transport and water level differen- 
ces by a step-by-step integration in time. In- 
stead of using such explicit equations, the 
implicit ones have been adopted; this has 
entailed gaining some advantage at the bound- 
ary of our region, the Gulf of Bothnia, which 
is somewhat complicated in shape. 
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2. Differential Equations 


Hydrodynamical equations can be written 
in the form: 


ou ian ae du 
a NE 
Ov ae dv 

2 


were x, y, 2 are the Cartesian coordinates, z 
being calculated as positive from the undisturb- 
ed surface of the sea upwards. 


£ is the elevation of the sea surface 

t is the time 

u, v, w are the velocity components 

y is the vertical eddy viscosity coefficient 
f = 2w sin y, the Coriolis parameter. 


It is assumed that the basin is not too com- 
plicated in shape, and that the water is incom- 
pressible and homogeneous. 

In these equations, the term of horizontal 
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viscosity and the non-linear terms are omitted; 
it appears that their effect is relatively minor. 

The above equations can be integrated over 
the depth from the bottom at — h to the surface 
at &, thus 


¢ 
dU ac du 
= a avon fa 


-h 


5 

oV ae Ov 

ae (+05 -JU+r [Tee 
5 


& 
where U= fu dz, and 
uh 


(a 
V = [vdz are the components of ho- 
-h 
rizontal water transport. 


u 
The stress components at bottom » (=) 
-h 


and» (5) are replaced byrU andrV respec- 
-h 


tively, where r= Le assuming the validity of 


ah? 
steady state of EKMAN (1905). On using Ek- 
man’s friction terms it is assumed that the 
depth is small. 


As is usual, the wind shear components 


du and er have been taken as tT, = 
Oz t dz £ 


= VHaWa and Ty=yvaws The wind compo- 
nents are u, and v,, w, the total wind speed, 
and y a constant. 

Following substitution of the shear values, 
the final equations are obtained, and these, 
together with the continuity equations inte- 
grated over the depth, constitute the complete 
set of equations necessary for the calculations. 
In addition, we have that ¢<h; thus we can 
write h instead of h + Z, except in derivatives. 
The equations are: 


QU ae 

en wala Td Gh 
oV ac 

Ot cl aed at Neh Re 
ac ae 

dt ax dy 


It is our aim to transform these equations, 
in order to make them suitable for calculations 
in a quadratic network. It is not advisable to 
take the available difference equations which 
are the simplest possible: it is easier to introduce 
the boundary conditions if they are made as 


symmetrical as possible. Thus we take: 


n u At n n+I 
CE UE Usp) 
n n+1 gh n en 
(Vie aed ee ra 
n ntT 
RE rg +T,At 


oh 
N n n+I S n _ pn 
An AG k+ı RE; 
TH a ee À 
+o eat a + TyAt 
At 
n+I n n n 
Ci a ee us 
n+I n+I n n 
U te Has + ee Frs 
n+I yn+I 
Fe Fam 


where Af is the length of the time step used, 
As is the mesh width, 


n gives the number of the time step 


calculated, 


j, k give the grid points which corre- 


spond to the coordinates x and y. 


These equations are implicit in the unknowns 
q P 

+ < n+ 
ty Ken and ¢"**, and thus must be cal- 


ik? Sik? 


culated by an iteration method. When the 
calculations are being started for a given time 
step, say n+ 1, the unknowns on the right side 
of the equations will have as substitutes other 
values which are suitable, here the known 
values for the time step n. The values thus 
obtained from the equations are used in the 
following step of iteration, and so on until the 


desired accuracy is attained. 


The first step in the calculations is that of 
feeding the initial values of U, V and € into 
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Fig. ı a, b. Water elevations in cm after 6" wind of 10 m/s from a southwesterly direction and the corre- 
sponding water transports. 


the machine. The boundary conditions must 
be taken into account with each step of iteration. 
They are as follows: 


(i) Near the coast, the current must be parallel 
to the coastal line. 

(ii) At the open ends, the water elevations 
must be given. 


3. The Stability of the Model 


When computations are being made by elec- 
tronic calculators, an examination must be 
made of the stability or instability of the differ- 
ence equations. If the method of Lax and 
RICHTMYER (1956) is used, it is easy to demon- 
strate that our equations, without Coriolis and 
frictions terms, are stable. Nevertheless, practi- 
cal calculations with the machine showed that 
the bottom roughness, the boundary conditions 
and the Coriolis force were causing instability 
and the friction decreasing the instability. To 
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make the equations stable, it was necessary to 
add a smoothing term to the right sides. For 
the first equation, the smoothing term was: 


J, R+1 j+1,k n 


7 i 


(Out Pier? D + U° ) 
a kJ" 


Similar expressions were used for the other two 
equations. The added Laplacian does not influ- 
ence the difference equations tending to the 
differential equations, when As and At tend to 
zero. In order to climinate the instability 
occasioned by bottom roughness, « had to be 
chosen at as high value as 0.25. 


4. Numerical Results 


To clarify the phenomena of current and 
water level in the Gulf of Bothnia, knowledge 
of which is still rather inadequate, two princi- 
pal cases were examined: that of a constant 
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Fig. ı c, d. Water elevations in cm after 12" wind of 10 m/s from a southwesterly direction and the corre- 
sponding water transports. 


wind blowing lengthwise over the bay, and 
that of a similar wind blowing transversely 
across it. For the calculations, the mesh width 
As of 44 km was chosen, and the time step At 
as 10 minutes. With these values and the depth 
h = 200 m, the Courant-Friedrich-Lewy crite- 
rion 


er Les 
Vgh 5 <V2 


for the stability of a simple wave equation is 
fulfilled. The wind-stress constant was chosen 
as y = 1.9-10 6 (cf. FISCHER, 1959). 


I 


In the first case, a constant wind was blowing 
from the southwest in the direction of the 
axis of the Gulf of Bothnia. As regards the 
boundary condition at the open end, the sea 
level was assumed to be normal (i.e. zero). 


The initial conditions were such that sea level 
elevation and current velocities vanished. Thus 
the wind started to blow simultaneously over 
the whole of the region. The water transport 
became strong in narrow sounds, and was on 
the whole more marked on the eastern coast 
and weaker on the western coast, as a result of 
the Coriolis force. (See Fig. 1 d, e, f.) 

At the same time, the water in most of the 
southwestern part displayed a slight subsidence, 
and then rose steadily in the northeast except 
in the Quark, where there was at first also a 
subsidence in water elevation. On the eastern 
coast, there was a noticeable piling up of the 
water masses due to the Coriolis force, as there 
was in the Quark, where the tilting of the 
water surface was very marked. The water 
transport to the northeast was first stopped in 
the north, and in the Bay of Bothnia after a 
period of six hours from the start of the wind. 
(See Fig. 1.) At the same time, the sea level had 
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Fig. ı e, f. Water elevations in cm after 3" wind of 10 m/s from a southwesterly direction and the corre- 
sponding water transports. 


the most irregular shape, tending later on to 
become smoother. A survey of the situation 
after 12 hours showed a clockwise circulation 
in the Sea of Bothnia. This circulation would 
however disappear later, and after 37 hours, 
when the currents had become steady, a north- 
easterly water transport was to be seen along 
the coast, and a counter transport in the centre 
of the sea. The sea surface is not a plane, but 
the general character is that the tilting is 
stronger in such areas, where the sea is shallower 
(P. WELANDER, 1957). 

The calculated magnitude of the difference 
in water level elevations between the extreme 
ends of the Gulf of Bothnia was checked 
against an actual case on 16 December 1953. 
On the roth December a weak southwesterly 
wind grew to a value of 9 m/s and lasted for a 
few days. The observed difference of 29 cm 
compared well with a calculated difference of 
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32 cm for a wind of 10 m/s; a calculated differ- 
ence of 26 cm would have been obtained if a 
wind velocity of 9 m/s had been used in the 
calculations. 

In the basins there occur regular oscillations, 
seiches. They show themselves most clearly in 
the Sea of Bothnia, where oscillations between 
the Swedish coast (Gävlebukten) and the 
Finnish coast can be observed. (See Fig. 2.) 
According to our calculations, the oscillations 
have a period of 5.45. This oscillation time was 
checked against that obtained from the formula 


I 
FES au 
J Veh 


of du Boys (Foret, 1894), which yielded 5.35. 
Here l is the distance between the two points 
of observations. 
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Fig. 2. Water elevations as a function of time in the 

southwesterly and northeasterly ends of the Sea of Both- 

nia, and in the southwesterly and northeasterly ends of 
the Bay of Bothnia. 


The water elevation in the Bay of Bothnia 
also shows slow variations similar to seiches, but 
they are heavily damped and seem to be caused 
only by the filling of the basin. 

It is still to be observed that in our calcula- 
tions no seiches which belong to the whole 
basin have been obtained. 


Ss. UUSITALO 


An analysis of some mareograms from mareo- 
graph stations on the Gulf of Bothnia yielded 
some oscillations which were well in accord 
with the assumption of seiches. Nevertheless, 
such a large number of local disturbances and 
parasitical oscillations confused the picture that 
no definite conclusions can at present be drawn. 


II 


The development of the case in which the 
wind blows across the Gulf of Bothnia generally 
followed the same lines as the foregoing. (See 
Fig. 3.) The maximum irregularity was here 
attained 3 hours after the start of the wind. 
At this stage the water was running out from 
the Gulf of Bothnia, due to the fact, that the 
strait to Baltic lies to the southeast of the south- 
westerly end of the Gulf of Bothnia. After the 
situation had become stable, there was a cir- 
culation of water, but even though the phe- 
nomenon was similar to the first case, it was 
not so clearcut. The water transport was over 
shallow areas with the wind direction, but 
hampered here by the direction of the coast 
line. 

In this case seiches would also occur. The 
oscillation time in the Bay of Bothnia was 
found to be 3.0—4.0" and in the Sea of Bothnia 
about 3.8". This is much less than the real 
seich-periods, and thus one may expect that 
the oscillations found may be spurious and 
caused by the finite difference approximation. 


5. The Influence of Density Variations and 
Friction 


When the water is not homogeneous in the 
whole basin, some complications occur. To 
illustrate what can happen, we assume that we 
have initially a two-layer system with a hori- 
zontal interface. (See Fig. 4.) The water be- 
neath this interface is denser than the upper 
strata. For the sake of simplicity, we assume 
that the basin runs east-west and the wind 
starts blowing from the west. The water near 
the surface is then transported to the east and 
piles up there. The increasing pressure therefore 
pushes down the layers underneath, and con- 
trariwise, due to the decreasing pressure at the 
western side, the water rises there. In approach- 
ing the position of equilibrium, the interface 
is much more tilted than the surface of the sea 
because 
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Fig. 3 a, b. Water elevations in cm after 3" wind of 10 m/s from a northwesterly direction and the corre- 
sponding water transports. 


where Ah and Ah’ are the differences in elevation 
of the surface of the sea and of the 
interface, respectively 

o and o’ are the densities of the lighter and 
denser water respectively. 

This condition serves wholly to disturb our 
friction calculations. Instead of taking the total 
depth of the sea for the depth in computing the 
friction, perhaps the depth of the pycnocline 
should have been taken. The situation is indeed 
far more complicated because circulation. is 
taking place not only in the upper layer, but 
also in the water masses below the pycnocline. 
When this reaches the sea surface, the situation 
becomes even more intricate. 

When the density distribution is more com- 
plicated than described here, it is clear that 
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there is no simple method of correcting our 
calculations. 


The bottom friction has been calculated 
according to the Ekman theory, which is valid 
for a special kind of circulation. It is an easy 
matter to determine the size of the errors intro- 
duced in special cases. Let us assume that the 
wind is beginning to blow over an undisturbed 
sea surface. Then for different stages of devel- 
opment the velocity-depth curve can assume the 
forms given in Figure 5. The velocity curve 
which corresponds to the Ekman theory is also 
plotted. The abscissae of the dots seen in the 
figures give the computed frictions for different 
cases. The crosses, again, show the correct fric- 
tion. The separation between the dot and the 
cross in each of the several sets represents the 
error. As seen in the figure, the friction comput- 
ed can even have a sign which is opposite to 
the correct one. 
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Fig. 3 c, d. Water elevations in cm after 28" wind of 10 m/s from a northwesterly direction and the corre- 
sponding water transports. 


Even though this might seem to invalidate 
the method of calculation, the effect on the 
results is after all minor. This is because in the 
deep sea the friction is rather small, and in 
shallow regions, where the friction is notice- 
able, the law of friction is in rather close 
agreement with the Ekman theory. 


6. Combination of the two Principal Solu- 
tions 


Assume that at the time ¢’, when U= V= 
€= o a constant wind starts to blow over the 
whole region, giving rise to the constant wind 


stress components T,, Ty. Then the solution of 


the equation system is 
U st Un (bat) 0, (tat) 
V=t,V,(t-t')+7,V,(t-t) 
C=ty bn (t—t') +t, Cy (t-2), 


where Ux, Vay x and Uy, Vy, Cy are the solu- 
tions of this equation system when 
the wind-stress components 7%, Ty 
are replaced by 1, o and o, ı re- 


spectively. 


If the wind changes but remains uniform 
over the whole region, then the increments of 
the wind stresses and he corresponding solu- 
tions may be calculated and added together. 
(See Fig. 6.) We obtain 


us I [ar U 
V2 > [ar 
ns 


U,(t~4)] 
V(t -1) + Atw V,(E-5)] 


U, (t=) +47 U 


Ten alt Fy) eee ee 


These equations are suitable for practical 
calculations. 
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Fig. 4. Influence of the pycnocline. 
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Fig. 5. Influence of the law of friction employed. 
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Fig. 6. The summation of constant wind stresses to 
obtain variable ones. 


From the theoretical point of view, it can 
be noted that these solutions tend to follow 
limiting values when Ar is approaching zero. 


felt 


U,.(t —t)dt + 
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This method was proposed by P. Welander in 
May 1957.We hope to test 1t on some meteoro- 
logically interesting cases in the near future. 
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Abstract 


Results of C14 determinations in surface water from the Pacific were in agreement with those 
reported by Rafter and Fergusson. However, abnormal C!4 concentrations seem to exist locally, 
for which no oceanographic explanation can be given. It seems premature to draw conclusions 
from existing determinations as to the rate of increase of C14 in surface ocean water resulting 
from the uptake of C14 produced in the atmosphere by atomic bombs. 

Samples from a constant depth of about 3,500 meters show a C14 content decreasing from 
south to north. This decrease may be attributed to radioactive decay of C14 during the time of 
migration. From this the northward component of the rate of water movement of about 0.06 
cm/sec can be calculated. The C!? determinations, for the purpose of correcting the C!* values 
for isotope fractionation effects, were found to be remarkably consistent, although made on 


reburned acetylene. 
Introduction 


The C14 content of the bicarbonate of ocean 
water is of twofold interest: First, it offers a 
valuable method for the study of isotope ex- 
change and transfer of carbon dioxide between 
the atmosphere and the ocean. Second, it can 
be used as a new tool in oceanography for 
the study of the movement of deep ocean 
water and the mixing of water masses. The 
rate of transfer of CO, from the atmosphere 
into the oceans was first discussed in the light 
of the C!* distribution by REVELLE and Suess 
(1957), and independently by Arnoıp and 
ANDERSON (1957). An order of magnitude of 
about ten years was derived by these authors 
for the residence time of CO, in the atmos- 
phere before it gets dissolved in the oceans. 
Subsequently several others, Craig (1957), 
FERGUSSON (1958), BOLIN and Erıksson (1950), 
attempted to establish a more precise figure 
for this residence time. All these considerations 
were based on the measurement of two em- 
pirical quantities: the apparent C14 age of 
ocean water, and the effect of industrial coal 


combustion on the specific activity of atmos- 
pheric CO. But the data up until now have 
been scarce, the observed effects o {small magni- 
tude, and there has been uncertainty about 
such things as the size of the so-called geo- 
chemical carbon reservoirs (such as the terres- 
trial biosphere). However these data are 
treated, the calculated residence time can 
scarcely be very accurate. But the addition of 
artificial radiocarbon to the atmosphere by 
bomb explosions during the years from 1954 
to 1958 will offer a much better opportunity 
to determine the accurate residence time and 
the geochemical pathway of CO, than the 
previous data afforded. A further discussion 
of this problem will therefore not be indulged 
in. Much more reliable answers will result 
from the combination of the data reported 
here with future observations during the 
coming decade. This paper will therefore be 
limited to the consideration of the use of 
radiocarbon as a means of following water 
movements. 

Part of our data were obtained from samples 
collected during the International Geophysical 
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Locations of sampling stations. Circles denote 
places where samples were taken by the Scripps Institu- 
tion. Open circles indicate surface water; solid circles are 
deep water samples or profiles. Crosses show locations 
from which samples were analyzed by the New Zealand 
Institute for Nuclear Science. 


Year for the study of secular variations (see 
Figure 1). During the next few years the C14 
activity of surface water will undoubtedly 
increase continuously, as a result of bomb- 
produced C!* in the atmosphere. Measure- 
ments made at this laboratory and at the Insti- 
tute for Nuclear Science in New Zealand will 
be a useful basis for future study of this prob- 
lem in the Pacific Area. 

The question of normalization presents a 
serious problem in the comparison of data 
from different laboratories. The oxalic acid 
standard distributed by the National Bureau of 
Standards is now used widely. In order to 
correct for isotope fractionation effects, results 
are normalized to a 6 C13 of —25 per mil 
relative to the Chicago PDB standard (Craıc, 
1957 a). The Lamont normalization, which we 
have adopted, uses this correction as well as a 
factor of 0.95 times the 6 C14 of the oxalic 
acid standard (BROECKER and WALTON, BROEC- 
KER and Orson, 1959).! 


1 It should be pointed out that the reference C14 in the 
oxalic acid standard is not corrected for any deviation of 
its ) C!3 from the chosen standard value of —25 per mil. 
The 6 C!3 of the oxalic acid standard is actually —18.65 + 
0.13 per mil relative to the PDB standard (Private com- 
munication of Toshiko K. Mayeda to J. R. Arnold). In 
case of deviation from this value, due to isotope fractiona- 
tion in the laboratory, the Ô C14 measured for the standard 
should be corrected accordingly. 
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Sampling and measurements 


On the Downwind Expedition of the Inter- 
national Geophysical Vear, six surface samples 
and ten samples from deep water werecollected 
and the carbon dioxide extracted and processed 
directly on shipboard. The samples for the 
profiles were collected with the same sampling 
equipment but stored in so-gallon steel 
drums and processed on shore shortly after 
the return of the ships. The Atlantic samples, 
obtained from Dr. John Lyman, were in the 
form of barium carbonate and were collected 
and processed in a different manner. 


The water samples were taken in a sampler 
of approximately so-gallon capacity construct- 
ed in the Special Developments Shop of the 
Scripps Institution of Oceanography under 
the direction of Mr. J. D. Frautschy. The 
sampler consisted of a stainless steel barrel 
about five and one-half feet long and fifteen 
inches in diameter, closed at each end by 
stainless steel flanges sealed with soft plastic 
gasket material. When closed, the end flanges 
were held together in place by a strong steel 
spring in the axis of the barrel. The flanges 
could be opened by a detachable “cocking” 
lever and held open about five inches from the 
ends of the barrel by a tripping mechanism 
which would be released bya messenger. When 
lifted inboard, the water sample was transferred 
as quickly as possible through a rubber hose 
into a “process tank” on deck built from a 
so-gallon glass-lined hot water heater, and 
filled with nitrogen gas before the transfer to 
avoid contact with the carbon dioxide of the 
air. The tank was fitted with a tightly bolted 
top plate with lead-in tubes for the addition 
of acid and the inflow and exit of nitrogen 
gas to flush out the liberated carbon dioxide. 
Two liters of concentrated HC] was added to 
each charge of the tank, and nitrogen was 
passed in through a coarse diffusion head at a 
rate of from one to four liters per minute, and 
finally led into an absorption bottle through a 
sintered-glass diffuser. The process tank was 
warmed with heating coils to a temperature 
of 40 to 50° C, under which conditions it was 
ascertained that all the CO, in the water was 
liberated and absorbed in a minimum of four 
hours. Twice this time or more was used in 
routine runs. 


The carbon dioxide in the nitrogen stream 
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was absorbed in one liter of concentrated 
ammonia solution saturated with barium 
chloride. We believe that absorption in am- 
monium hydroxide solution is much faster 
than in alkali hydroxide, and that the use of 
ammonium hydroxide represents a real im- 
provement over techniques used elsewhere. 
The one-liter bottle, containing the remaining 
solution and the precipitated barium carbonate, 
was tightly stoppered and returned to the 
laboratory. This method proved very satis- 
factory, in spite of the fact that the precipita- 
tion of solid BaCO, on the diffuser head 
generally made its replacement necessary at 
least once during a run. Completion of the 
absorption could always be verified, however, 
by replacing the absorption bottle with a 
fresh one and observing whether precipitation 
of barium carbonate continued or not. In the 
absorption of the CO, from some of the later 
samples, strontium chloride was substituted 
for barium chloride. In any case, the strontium 
or barium carbonate was reduced with magne- 
sium metal to the carbide, which was then 
hydrolyzed to acetylene. 

Hydrochloric rather than sulfuric acid was 
used for acidification with the thought that it 
would interfere less with the subsequent deter- 
mination of radium. In five liters ofeach sample, 
after the removal of carbon dioxide, barium 
sulfate was precipitated by the addition of 
barium chloride. This precipitate, carrying 
with it the radium, was washed and separated 
by decantation and returned for the determina- 
tion of radium. 

The deep samples for the C14 measurements 
were all taken from depths of about 3,500 
meters. In each case a Nansen bottle was 
suspended about five meters below the sampler 
in such a way that the bottle, with its protected 
and unprotected thermometers, was tripped 
when the large sampler closed. From these 
thermometer readings the depth was accurately 
determined in the standard manner. 

With one exception, the surface samples 
were obtained through the ship’s salt-water 
pumping system, which operates without 

iffusing air into the water. One water sample 
was obtained from mid-depth, 400 meters. 

The La Jolla Radiocarbon Laboratory uses 
the method developed at the U.S. Geological 
Survey in Washington (Suzss, 1954) in which 
acetylene serves as the counting gas. The C14 
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counting was done in a special low-level counter 
developed by HourerMans and OESCHGER 
(1958) and manufactured in the instrument 
shop of the Physikalisches Institut, University 
of Bern, Switzerland. This counter showed a 
background of less than two counts per minute 
at one atmosphere and a sensitive volume of 
about 1.7 liters. Total counting time for each 
sample was several days. 


Discussion of results 


The results of the measurements are given 
in Table 1, 2 and 3, and shown in Figure 2. 
They contain two remarkable features: (1) 
the occurrence of abnormal surface water 
concentrations of C!4; and (2) a distinct, 
regular increase in the apparent age of deep 
water from south to north. 


Surface water samples. Careful measurements 
in the laboratory of the Institute for Nuclear 
Science, New Zealand (RAFTER and FERGUSSON, 
1957; FERGUSSON, 1958) have shown that the 
activity of the surface water in the South Pa- 
cific decreases from lower toward higher lati- 
tudes along with decreasing stability of the 
thermocline. Values more than 150 per mil 
below the equilibrium value are found beyond 
the Antarctic convergence. Considering the 
uncertainty in the calibration and the experi- 
mental errors of about one per cent, our results 
agree very satisfactorily with the New Zealand 
findings. In comparing results we have assumed 
a difference between the New Zealand and 
Lamont normalizations of -61 per mil. In a 
few cases, however, there is considerable 
deviation from the general trends. One ab- 
normal sample was found by Rafter and Fer- 
gusson (BURLING and GARNER, 1959) with an 
activity of more than three per cent below 
the expected, and the investigators attribute 
this deviation to chemical effects from storage 
in a steel drum. We believe that in our case 
such an explanation of the deviations is highly 
improbable. Along the La Jolla coast and in 
offshore areas varying AC14 values ranging 
from — 50 to about — 80 per mil have been 
obtained. These variations can be explained 
by upwelling water along the California coast. 
It is well known that such upwelling results 
in great differences in the local properties 
of surface water at different times and locations. 
The vertical circulation and mixing of surface 
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Fig. 2. Upper part: AC! plotted as a function of latitude for surface water samples (open 
circles indicate La Jolla measurements, crosses New Zealand measurements), and for samples 
from approximately 3,500 m. depth (solid circles). 

Middle part: Normalization and the conversion of AC!4 to apparent ages is indicated. All 
values, except that for the NBS oxalic acid standard, are normalized to a ÔC13 of —25°/00 
relative to the Chicago PDB standard. 

Lower part: AC!4 as a function of depth. Crosses indicate New Zealand measurements. Open 
circles show results on samples collected by J. A. Knauss, and solid circles by J. E. Reid. 
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Table 1. IGY downwind expedition’ ~ 


Collector: N. W. Rakestraw 
Dates of Collection: October through December 1958 


C8 measurements: P. Eberhardt in laboratory of H. C. Urey 
Seen ee RE ee = 


: 2 ol ue 
o ° = S te 5 
: = = Bost le 
> © = = > nm SL dan 
5 a 3 2 2 > 2) Baa EEE: 
eee 3 & | £ |=] o [gas are 
| Si S =) Ö e) = Sa Lee 
58 2 1)” SRY IN 125, W o 34.69 4.82 24.4 —0.5 —40 
62 6 5° 7308 W o 34.40 4.65 28.0 —0.1 —40 
66 8 TRES 1325 W o 34.68 4.88 26.7 —1.0 —42 
68 II SAMEORS LS SES NI fo) 35.05 es ET; —0.4 a 
go I 46° 30'S 1103 W o 7 9.2 —0.7 —60 
93 be 125 ES 96° W o 6.64 11.5 —0.6 —62 
60 4 Tr NE 0728737 400 34.65 0.13 10.9 —1.5 —103 
59 3 14° 58’ N TITTEN 2987 34.69 2.73 1.66 —I.2 —230 
61 5 708 «IN| 129° 16’ W 3508 34.69 3.43 1.48 — 1.4 —243 
63 7 2 085 | 131 27 3473 34.69 | 3.48 1.52 15 —235 
57 9 1 9368 135° 30 3449 34.69 3.73 1.58 —1.8 —220 
67 10 25° S 145° 3560 34.65 4.05 1.54 —2.9 —204 
69 12 1023548 1324 10 3500 34.67 4.00 1.50 —1.2 — 195 
88 13 46° 46° S 1290053 3550 34.69 4.03 1.62 —2.0 —200 
OI 15 AO TES 1030.20 3500 34.67 4.10 1.85 —2.5 —180 
94 m7 SOMMES AS 84° 10 3500 4.15 1.46 19 —196 


1 The values tor AC given here differ slightly from those previously reported (Rakestraw, Oeschger 
and Suess (1959) Preprints, International Oceanographic Congress, p. 440). They are calculated 


with a revised calibration. 


Table 2. Samples from La Jolla beach or offshore 


13 0 14 0/ 
LJ No. Material | Le = Behe Date of Collection 
86 Shells, collected alive —8.3 —85 1953 
(mytilus Californianus), 
La Jolla Shores 
48 Vellela velella, SW of Farallon —23 —82 15 April 1957 
Islands 
70 Kelp (Macrocystis pyrifern) —16 —70 8 July 1957 
65 Red Tide (Gonyaulax) —18 —60 18 August 1958 
97 Shells, collected alive + 0.6 —37 February 1959 
(mytilus Californianus), 
La Jolla Shores i 
127 Carbonate from Pier Water — 2.3 —88 4 August 1959 


water with that from lower depths results 
in variations in the C!* activity of surface 
water as well as in other characteristics, such as 
temperature. One measurement (LJ 149) gave 
an abnormally high value for which we have 
no explanation. 

As Fercusson’s results show (1959) one 
must be cautious in interpreting increases in 
the C!* activity of surface water as a conse- 


quence of uptake of bomb C14 from the 
atmosphere. Many measurements will be 
necessary to determine the rate of this uptake. 

Deep water samples. The C14 concentration 
of a particular deep water sample may be in- 
fluenced by the following factors: 

1) The degree of equilibration between the 
CO, in the water and that in the atmosphere, 
provided the water mass can be regarded as a 
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Table 3. Profiles 
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closed system. If not, then the water mass 
must be considered as a sum of components, 
each with a different history. 

2) The time since such equilibration occurred, 
during which radioactive decay has taken 
place. 

3) Oxidation of particulate or dissolved 
organic matter. Settling of particulate organic 
matter from the ocean surface and subsequent 
oxidation will add younger carbon and 
decrease the apparent radiocarbon age. Evi- 
dence of sucha process may possibly be obtained 
from the O, distribution in the water. 

4) Redissolving of carbonaceous sediments 
from the ocean bottom. Geological evidence 
indicates that such processes are in general 
too slow to affect C14 concentrations. 

The Pacific deep water obtained its C14 
through exchange with the atmosphere south 
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of the Antarctic convergence below 60° S. 
latitude, or through admixture of surface water 
from lower latitudes. New Zealand measure- 
ments show that exchange is incomplete and 
the water masses sink with a AC! of about 
— 150 or less (apparent C14 age of 1,200 to 
1,500 years). Our measurements show a further 
increase in apparent age of 200 to 300 years 
as the water moves toward the equator. This 
age difference may well indicate the true 
travel time of the water mass. At the same 
time a decrease in O, occurs that can be 
attributed to oxidation of organic matter. 
Assuming that this organic matter reached 
the deep water by sedimentation from the 
surface and is correspondingly high in C14, 
this addition of younger carbon to the dissolved 
bicarbonate will make the radiocarbon age of 
the water appear too young. Sea water con- 
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tains about 2.3 millimol per liter of carbon 
in the form of bicarbonate. On its way from 
40°S to 15°N the water loses about 1.3 
milliliters or 0.055 millimol of oxygen, equiv- 
alent to about 4.2 % of the carbon present. 
The carbon added by oxidation of organic 
matter may be about 20 per cent higher in 
C14 and therefore may increase the C!* con- 
centration of the water mass by 0.042 x 20 % 
= 0.84 %. This corresponds to a decrease in 
the apparent age of about 100 years. Figure 2 
shows that the difference in the apparent C14 
age for the samples from around 3,500 meters 
depth between 40°S and 15°N was found 
to be between 200 and 300 years. If the 
‘oxygen in this water mass was consumed by 
organic carbon from the surface, then the true 
travel time of the water would be about 100 
years longer. Assuming that effects from other 
processes, such as the dissolving of particulate 
inorganic carbonate, are unimportant, then 
we can estimate the travel time for 55° of 
latitude (6,100 km) as 300 + 100 years. This 
corresponds to an average northward com- 
‘ponent of the velocity of the water of 0.06 + 
0.02 cm/sec. 


. Our evidence leads us to expect that radio- 
‘carbon studies will be more valuable in ex- 
‘plaining water movements in the Pacific than 
ain the Atlantic, where conditions are appar- 
‘ently more complicated. 


Profiles. The variation of the C!* concen- 
tration as a function of depth in the Pacific was 
investigated at two stations. The samples from 
one station were collected by Dr. John Knauss 
on the Doldrums Expedition of the Scripps 
Institution of Oceanography; those of the 
other station by Mr. Joseph Reid on another 
cruise. In both cases the seawater samples were 
brought ashore in 55-gallon drums and proc- 
essed at La Jolla. The results obtained show a 
more regular depth dependence than those 
obtained by the New Zealand laboratory. 
The C?!4 concentration was found to decrease 
almost linearly with depth down to about 1,000 
meters and from there downwards to remain 
almost constant (Figure 2). One profile 
showed evidence for a C!4 minimum at a 
depth of about 2,000 meters. Until the occur- 
rence of this at other stations is verified, specula- 
tion as to its significance is scarcely profitable. 
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C13 results. The mass spectroscopic C1? deter- 
minations carried out by Dr. Samuel Epstein 
and Mrs. Iren Vidzuinas at the California 
Institute of Technology, and in Dr. Urey’s 
laboratory at La Jolla by Dr. Peter Eberhardt 
and Dr. Robert Bowen, are remarkably con- 
sistent. They show a relatively small spread, 
much smaller than that found at the Lamont 
laboratory (BROECKER, TUCEK, and OLson, 
1959) in samples from the Atlantic. Also, their 
variations follow qualitatively the expectation 
that the C13 should show lower values at 
depths where oxidation of organic matter 
may have taken place, than at the surface 
where photosynthesis predominates. The lowest 
C18 value observed (LJ 153) is in a sample 
which comes from a depth near that of the 
oxygen minimum. However, it should be 
noted that the C1 determinations were made 
on reburned acetylene, and that small isotope 
fractionations in the laboratory during extrac- 
tion of CO, and preparation of acetylene 
cannot be excluded. 
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On the Frequency of Lightning Flashes to High Objects 
A Study on the Gulf of Bothnia 


By D. MULLER-HILLEBRAND, Uppsala 


(Manuscript received May 4, 1960) 


Abstract 


The space-charge calculation of a lightning stroke according to Golde results in a simple geomet- 
rical relation between the area within which a flash strikes an object and its height. A comparison 
between the number of lightning strokes to chimneys and the number of lightning downstrokes per 
km? shows, that the statistical attraction distance of a 70 m high chimney is more than 250 m 
and thereby more than the leader’s step-length after Schonland. An explanation may be the 
space-charge plume, which leaves the lightning rod in a thunderstormfield and may conduct a 
lightning discharge to the chimney from some hundreds of meters’ distance. 


High objects and leader head” 


In 1933 B. WALTER in a report which 
attracted much attention at the time pointed 
out that the lightning flash in its approach to 
the earth’s surface does not fix on its down- 
stroke point until the final phase at a distance 
of less than 150 m from the ground. Walter’s 
statistical work was based on 438 cases of 
lightning damage in Hamburg during 27 
years. The average number of cases of light- 
ning damage annually is calculated at 0.57 
per km? in Hamburg with ca 19 thunderstorm 
days a year. The area A, within which the 
flash strikes an object with a height h above 
its surroundings, is calculated from geometrical 
ratios as 


Fe (= a 1) (th) Mens 
A = 27h" d (d>h) (x a) 
A = den (d<h) (1b) 


with d as attraction distance. The area and in 
consequence the probability of attracting a light- 
ning stroke increases linearly with the height 
of a lower object and is independent of the 


height of full objects. This simple relation has 
been the basis of many investigations. R. H. 
GOLDE (1945) has presented a plausible opinion 
that the attraction distance is dependent on the 
size of the electric charge in a stepwise pro- 
jecting space-charge channel, which prepares 
the lightning path. The distance is greatest 
for large electric space-charges and consequent- 
ly also for high lightning current-strengths 
(Schonlandprocess, step process of the initial 
leader). 

Golde’s pre-requisite condition for the cal- 
culation is the assumption that the space-charges 
are distributed in a kind of “pipe” with a 
charge density which decreases exponentially 
with the height. S. B. Griscom (1958) puts 
forward views which lead to the assumption 
that the charges at the “leaderhead” are 
chiefly hemispherically distributed. The form 
of the corona envelope earthward of the leader 
is like an onion. This leaderhead contains a 
very large concentration of immediately avail- 
able energy. The leader grows by fits and 
starts in the direction of the ground, when a 
transformation takes place from the hemisphere 
to a channel and from the channel to a hemi- 
sphere. The distance from which a flash-over 
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takes place between the leaderhead and a 
high object is assumed to be 50 m. According 
to R. H. Golde, the smallest distance from 
which the flash starts its last step to the ground 
or from which a discharge from the ground 
to the space-charge of the flash starts is less 
than 17 m. Practical observations have con- 
firmed this. At a distance of less than 15 m 
from the ground the flash has altered its path 
and struck an object from the side (4, 5). In 
calculations of the area of attraction a distance 
of 40—60 metres from the ground is often 
assumed, 


Lightning flashes to chimneys 


For a period of 17 years from 1943 to 1959 
we have had under observation about 37 
chimneys at 27 different places. Twenty chim- 
neys (height 60—106 m) are situated on the 
coast of the Gulf of Bothnia (Fig. 1); 15 north 
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Fig. 1. Lightning strokes to 16 chimneys on the coast 
of the Gulf of Bothnia and number of lightning down- 
strokes per 100 km?. 
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of Gävle (60°40° lat., 17°02’ long.) to Hapa- 
randa (65°50° lat., 24°09’ long.) and 4 in the 
Gävle district. One was situated to the south 
at Karlskrona (s6°10’ lat., 15°30’ long.). The 
remaining 17 chimneys are situated inland. 
The observation was carried out by means 
of steel rods. Through the courtesy and as- 
sistance, continued for many years, of the 
managements of the respective factories, we 
have obtained information on the annual 
mounting of the steel rods from inspection 
cards. The average number of observation 
years is more than 12 and the number of 
measured flashes with more than 5 kA current 
is 28. The risk, of chimneys of an average 
height of ca 70 m being struck by a flash with 
a current strength of more than 5000 A, is 
3.83.90 (see Table I). Theis an average 
value for all observation years together (7.4 % 
for the years 1943—1949, richer in thunder- 
storms, and 4.7 % for the years 1950— 1958, 
poorer in thunderstorms). In Fig. 2 the results 
are reproduced graphically. The probability 
of the downstroke of a flash having a current 
strength of 20 kA and more is ca 2.7%. A 
powerful flash with a current strength of 50 


Table I. Comparison of the measurements on 
chimneys. 
Chimneys 
situated 
un inland 
coast 
Average number of thunderstorm 
daysapemwsyeatnL ee 5 12 
Number OL Chimneys rs ze ren 20 18 
Nwerasesheicht em: sameeren tae 70 60 
Number of observation years..... 1240001270 
Number of flashes measured 
oe enr dr re moe 18 16 
Over ASCUTTENt Er. ee 13 15 
Strike probability, % 
For flashes with at least 5 kA 
CUITE ee Ce 5.25 6.5 
Area of attraction, km? 
Path of the flash influenced with 
a distance Of = 100m 7. 0.026| 0.024 
CRD ON ES 0.084| 0.071 
Lightning downstroke per km? per 
year 
| Rope ce itCYON tad ee. POP || BGS 
(fh Sat PLOY NV aan Chen ame 0.625] 0.915 
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Fig. 2. Probability of currents with lightning strokes to 
chimneys. 


A Mean value 
C On shore 


B In the country 
D According to Hylthén-Cavallius 
and Strömberg (1959) 


kA and more will strike a 70-metre-high 
chimney annually with ca 0.75 % probability. 


If the area of attraction is calculated from 
formula (1) with the proportionately high 
value for the attraction distance d = 100 m, 
the area of attraction is therefore 0.0263 km? 
large for a chimney with 60 m effective height. 
Here it has consequently been assumed that 
the lightning path is influenced at a height 
of 60 m and at a distance of 100 m from the tip 
of the lightning leader. It was assumed in the 
calculation that the effective height of the 
chimney was reduced by 10 m by surrounding 
buildings. 

With A= 0.0263 and with a probability 
0.0525 (Table I) the average value of the 
annual number of lightning downstrokes per 
0.0525 
0.0263 
sponding value inland is 2.75. This is not 
likely. The number of years of observation 
with lightning counters is certainly not suffi- 
cient for drawing definitive conclusions (6). 
However, the following can be asserted. 


2 = 1.99 at the coast. The corre- 
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In 1959 the number of thunderstorm days 
on the east coast (Karlshamn to Haparanda) 
was 7 to 3. The average value for 10 years is 7 
in the south and ca 4 in the north, and the 
average value for 37 years is 8 and 5 respec- 
tively. During 1959 the number of lightning 
days correspond to a medium value. In 1959 
there was recorded an average value of 38 
lightning strokes to ground per 100 km? with 
14 lightning counters placed on the east coast 
(Karlshamn to Luleä). The highest number 
of flashes for 100 km? was 75 (the Stockholm 
archipelago) and the smallest number was 7 
(Umeä). If we count on an average annual 
value of 60 flashes per 100 km? on the east 
coast, then this will presumably be a high 
amount. Thus the attraction distance for 70- 
metre-high chimneys at the coast will be 309 m. 

The chimneys inland are situated north of a 
line from Gothenburg (57°40’ lat. 11°55’ long.) 
to Norrköping (58°35 lat. 16°10’ long.) in 
areas with 8-12 thunderstorm days. as an 
average value for 10 years and 8-13 thunder- 
storm days as an average value for 37 years. 
The lightning counters set up in the neigh- 
bourhood of the factories recorded in 1959 50 
flashes to the ground for 100 km?. The highest 
number was 110 flashes per 100 km?. This is 
comparable with the average value of 0.915 
reproduced in Table I, when an attraction 
distance of 250 m is used as a basis. 

The values reproduced in curve D in Fig. 2 
have been taken from a nine-year investiga- 
tion (1950 to 1958) by N. HYLTHÉN-CAVALLIUS 
and A. STRÔMBERG (1959), which comprises 
1125 observation years and 19 current values 
more than 5 kA. In the transfer of their results 
we calculated with 75 % of the number of 
station years (845), partly with regard to 
a 6% loss of measuring rods and partly 
(after consultation with the authors) with 
regard to the possibility that not all magnetic 
rods were mounted during the whole thun- 
derstorm season. After the reduction the 
probability values agree for flash currents 
over 35 kA (curve D and curve A). On the 
other hand, the probability of a flash down- 
stroke with a current strength of at least 5 kA 
is 2.75 % according to curve D (5.7 % curve 
A). The average height of chimneys struck 
by a lightning was 55 m and the effective 
height probably 45 m. The medium value of 
lightning strokes to earth is estimated at 


Tellus XII (1960), 4 


ON THE FREQUENCY OF LIGHTNING FLASHES TO HIGH © BPE CaAgs 


0.5/km? in the region with the majority of 
chimneys. The attraction area becomes ac- 


à 0.0275 
cording] ae 0.055 km? and the attrac- 
tion distance 212 m with h = 45 m. 

The result of the simple calculation is: the 
distance of attraction is about 200 to 300 m 
and consequently not in agreement with the 


Schonlandprocess in the interpretation of 
Golde. 


Corona discharges and space-charge 


For the regions of the U.S.A. rich in thun- 
derstorms (30—40 thunderstorm days per 
year) McCann (1944) has reported the strike 
probability for fire-watch towers (h,, = 28 m), 
for high radio masts (h„=107 m) and for 
chimneys and high buildings (171 m). These 
values and the ones measured on the Empire 
State Building in New York have been re- 
duced to 10 thunderstorm days (for lack of 
knowledge of lightning downstrokes for one 
km?) and reproduced in Fig. 3, together with 
our measurements. The probability increases 
approximately with the square of the height of 
the object. 
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Fig. 3. Strike probability for towers and chimneys. 
1. Empire State Building, New York 
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Fig. 4. Electric field strength and corona current on a 
7o-m-high chimney and a 15-m-high point. 


The higher a chimney is, the greater is the 
corona current which is generated from the 
lightning rod interceptor in the electrical 
thunderstorm field of some thousands of volts 
per metre with a thunderstorm. Fig. 4 repro- 
duces as an example the electrical field strength 
which was measured in Uppsala on 25 August, 
1959, in a widespread cumulo-nimbus without 
discharges. Some twenty metres from the 
field-measuring instrument there was a I5- 
metre-high point, the corona current of which 
was simultaneously recorded. At a distance 
of 5.85 km there was measured the current 
through the lightning rod on an 70-metre- 
high chimney. The current generates a posi- 
tive space charge at positive electric field 
strength. The current strength has a peak 
value of 250 wA; the charge generated for 10 
minutes is $2 m Coul. A positive current of 
250 uA leaves the rod and flows upwards or a 
negative current from the direction of the 
cloud to the interceptor flows to the ground. 


R. Davis and W. G. STANDRING (1947) 
have measured the discharge currents in the 
flying cable of a balloon, flown at 600 m height. 
In thundery weather they got currents of 
about 25 mA max. value during several min- 
utes. The recorded currents indicate the 
passage to earth of a large quantity of charge, 
3 to 4.5 Coul. The voltage gradient was also 
measured, but on the approach of a charged 
cloud difficulties occurred and no records 
could be obtained. The order of magnitude 
of the potential gradient may have been the 
same as with our measurements, and we can 
see that currents and charges increase approxi- 
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Fig. 5. Lightning stroke to a 45-m-high chimney (taken 
by R. Kinnunen at Thomajärvi in Finland). 


\ 2 
mately with (height)?, equal to (=) =74 


times the current of 250 uA=18.5 mA, anda 
charge of 52 m Coul. =3.85 Coul. Theoretically 
it is well known that the field strength of a 
point in a homogeneous electric field (long 
ellipsoid homogeneously polarized) increases 
with the square of the height. We get the con- 
clusion that the space-charges of these high 
objects affect the probability of a lightning 
stroke and that the probability increases with 
the square and not linearly with the height 
according to McCann (1944). 


Fig. 5 reproduces a lightning stroke to a 
45-m-high chimney, taken by R. Kinnunen 
at Thomajärvi on 5.7. 1954 at 8—9 p.m. 
(62°13° lat., 30°15’ long.). The lightning 
channel goes obliquely over a length of at 
least 145 m, possibly essentially more. The 
following explanation seems to be plausible. 
A space-charge plume leaves the chimney. 
The space-charges are spread rapidly due to 
their own electric field and as a consequence 
of the eddydiffusion in the atmosphere. It is 
well known that the dispersal of pollution 
from chimneys depends in a high degree on 
the micrometeorological conditions. With a 
low turbulence the smoke of a chimney drifts 
away as a narrow cone which does not diffuse 
until it has reached a considerable distance 
from the site (10). If a large number of con- 
densation nuclei or fogdroplets is present, the 
mobility of the ions decreases, and the electric 
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field strength on the edges of the space-charge 
may achieve values of a hundred kV/m. This 
changes the original electric field drastically: a 
lightning discharge may be conducted to the 
“plume” and thereby lead discharges from 
some hundreds of metres’ distance to the 
chimney. Presumptions for this are seldom: 
the probability of a stroke for a 70-m-high 
chimney may be about one in 100 lightning 
days on the coast and one in 130 to 150 
lightning days inland. 

The pre-requisite conditions are obviously 
not always present and are more easily ful- 
filled the higher the chimney. The existence 
of this positive space-charge at negative field 
strength is evident when the primary electric 
field between a thunder-cloud and earth 
disappears through a lightning discharge 
between the clouds, and the space-charge 
field alone is effective. A “rebound” discharge 
between the positive space-charge cloud and 
the lightning conductor may occur. Positive 
currents of 0.5 to 1.5 kA in magnitude were 
recorded by N. HYLTHÉN-CAVALLIUS and À. 
STROMBERG (1959) approximately as often as 
the number of lightning currents over 5 kA. 
Similarly McCann (1944) has several times 
measured currents of the same order of magni- 
tude (100-800 A, 50-150 usec duration) 
and their origin is explained as induced currents. 
However, it seems to be possible that these 
positive currents arise from the positive-charge 
cloud generated from the high object itself. 


An argument for this opinion in Figs. 6 and 


Fig. 6. Lightning flashes, starting from both towers, 

with upwards directed ramification. Tower distance 

400 m, in direction of the figure 260 m. Visible channel- 
length 850 m. 
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Fig. 7. Oblique discharge 6 min. 20 sec. later. The dotted 
line is the contour of the mountain, the top of which is 
concealed behind a cloud. Dischargelength at least 1.000 m. 


7, lightning flashes on the 70-m-high towers 
on Monte San Salvatore near Lugano, Switzer- 
land. We get these pictures in connection 
with lightning researches in collaboration 
with Prof. Berger, Ziirich. The exposure-time 
of a 16-mm-film was 8 sec. The discharges 
start from the top of the tower and the 


channel’s ramification is directed upwards. 


Experimentally this phenomenon is not 
confirmed owing to the difficulties in getting 
sufficiently large space-charges. But with nega- 
tive space-charges produced by impulse dis- 
charges we can get a back-discharge from 
these space-charges to the electrode. As an 
example Fig. 8 shows a negative corona from 
a 6.2 cm sphere with a 1/50-usec wave and a 
voltage of about 1 MV in a photograph kindly 
offered by Mr. O. Salka. From the space-charge 


flashes a spark goes back to the electrode. 
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Fig. 8. Negative impulse corona with a back-spark 
(photograph by O. Salka). 


This experiment is not a proofofthe existence 
of this phenomenon and may only give an 
explanation and an idea of how we imagine 
the occurrence of small positive lightning 
currents of the order of some hundreds and 
thousands of amperes with chimneys. 
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Phase Changes in the Daily 


Variation of the Cosmic Ray Nucleonic Component 


By N. R. PARSONS, Fysiska Institutionen, Uppsala 


(Manuscript received February 22, 1960) 


Abstract 


Results of the analysis of solar daily variations in the cosmic ray nucleonic component at nine 
stations during the International Geophysical Year are presented. 

Examination of month to month phase behaviour leads to the conclusion that there is an im- 
portant component variation synchronized in Universal Time at all stations, and its origin is 
discussed. Its effects on inferences about the primary anisotropy responsible for the main part of 
the daily variation and on various characteristics of the observed variation are also discussed. 


I. Introduction. 


It is well known that the average time of 
maximum. of the solar daily variation in cosmic 
ray intensity, measured by recorders at ground 
level, undergoes long term changes. The 
results of THAMBYAHPILLAI and ELLIOT (1953), 
SARABHAI et al. (1954), VENKATESAN and DATT- 
NER (1959), and others, have shown that these 
long term (e.g. year to year) changes are 
associated with the cycles of solar and geo- 
magnetic activity, although these associations 
are not entirely consistent and do not give 
a clear indication of the mechanisms producing 
and directly controlling the daily variation. 

Its local time dependence, together with 
other features, suggest strongly that the major 
part of the daily variation arises from an 
anisotropic distribution of the primary radia- 
tion entering the geomagnetic field region. 
This interpretation is now fairly generally 
accepted, as is also an explanation of the long 
term phase shifts in terms of changes in the 
mean direction of the primary anisotropy. 

A recent discussion of probable and possible 
causes of such an anisotropy has been given 
by DATTNER and VENKATESAN (1959) and it 


seems clear that several separate effects may 
contribute, subject to more or less independent 
modulation by changing conditions in inter- 
planetary regions. We shall not be concerned 
here directly with this problem but rather 
with an examination of the degree of agree- 
ment between different recording stations on 
the direction of the effective anisotropy and 
on changes in this direction inferred from 
observed phase changes in the daily intensity 
variation. 

Several authors have examined records of 
the meson component and shown that there 
is a high degree of correlation between the 
phase trends in yearly mean diurnal variations 
at different stations (cf. VENKATESAN and DATT- 
NER, 1959). With reasonable assumptions 
about the energy dependence of the aniso- 
tropy, sufficiently accurate allowance can be 
made for deflections of primaries in the 
geomagnetic field, to produce fair agreement 
on the mean direction of the anisotropy in 
different years (cf. DORMAN, 1957). However, 
the agreement between stations concerning 
shorter term (e.g. month to month) phase 
changes is not as satisfactory. In many cases 
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the apparent inconsistencies are considerably 
greater than could be expected from the 
relatively larger statistical errors inherent in 
the shorter term measurements. A similar 
situation exists for the amplitude of the ob- 
served daily variation. 

In the case of the meson component to 
which most of the published results refer it is 
tempting to attribute much of the incon- 
sistency to important residual atmospheric 
contributions. Adequate correction for the 
influence of changing temperature-height dis- 
tributions in the atmosphere cannot be made 
at present because of our very limited knowl- 
edge of these changes throughout each day. 
In periods as short as one month these unknown 
residual effects could add an important com- 
ponent to the daily variation. Such a compo- 
nent is certain to differ considerably between 
stations in any one month and to change from 
month to month, thus seriously affecting 
phase and amplitude correlations. 

Such influences however may not account 
fully for the inconsistencies. By examining 
records of the nucleon instead of the meson 
component these difficulties can be almost 
completely overcome since in this case the 
only significant atmospheric effect is due to 
total pressure changes for which correction 
can be accurately made. There is a further 
advantage in that the amplitude is in general 
larger than for the meson component and 
the accuracy of phase determinations is corre- 
spondingly greater. As will be evident from 
results for the nucleonic component presented 
in a later section, month to month phase 
changes at different stations remain rather 
poorly correlated, indicating that disturbing 
effects of non-atmospheric origin must also play 
a significant role. 

It is customary to resolve the observed 
daily variation into its harmonic components. 
The first two (diurnal and semidiurnal) usually 
suffice to give a reasonably close fit to the 
observed variation although it cannot be in- 
ferred that these have independent physical 
significance. For an anisotropy of any partic- 
ular form the actual course of the variation 
will depend on the directional acceptance 
pattern projected outside the geomagnetic 
field for each recorder, and several harmonic 
components will in general be required to 
describe this course accurately. The relative 
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magnitude of the harmonics will also differ 
from station to station. However, in the great 
majority of cases the diurnal component is 
found to be the largest and must reflect the 
existence of a single main “core” of the 
anisotropy. Its time of maximum will then 
in general give a good indication of the direc- 
tion in asymptotic longitude of this “core”. 
It has often been demonstrated (e.g. VENKA- 
TESAN and DATTNER, 1959) that the diurnal 
component does in fact exhibit much more 
consistent behaviour than the smaller semi- 
diurnal component and we shall confine our 
attention mainly to it. 

In what follows we shall present results of 
analysis of the daily variation in nucleonic 
component intensity at nine stations covering 
the period of the International Geophysical 
Year (July 1957—December 1958). A proce- 
dure will be described by which a significant 
general improvement can be made to the 
phase agreement between stations on a month 
to month basis. The hypothesis is made that the 
observed daily variation contains a diurnal com- 
ponent arising from synchronized Universal 
Time variations effecting all stations. In separate 
monthly periods such a U.T. component is 
determined (by relatively crude methods) 
which on subtraction from all stations yields 
the best possible improvement in phase 
agreement. The relative success of this proce- 
dure suggests that the U.T. component is 
genuine, and its effects on inferences from the 
observed daily variation will be briefly dis- 
cussed. 


2. Analysis of nucleonic component records 


Choice of stations and selection of data 


For any comparison of the mean daily 
variation at different stations, particularly when 
averages over periods as short as one month 
are considered, it is very desirable that exactly 
the same days should be used for all stations. 
Considerable changes in the form of the 
variation occur from day to day and inclusion 
of some days for which data are not available 
for all stations can only detract from the value 
of the comparison. On the other hand if a 
comparison is attempted between too many 
stations, one finds relatively few days with 
complete records available from all stations. 
Thus in the present case the analysis is limited 
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Table 1. Stations used in the analyses. ~ 


Dene ee 


Geographic Geomagnetic A ten 
Station 
Lat. Long. Lat. Long. (m) morte 
RU sp ee i ee u e 
lalelonee ten. SIP? BEM PA oma GAS | 555 225 725 7/57—12/58 
Uppsala re" tact: +59° 51’ 17° 55 E | +58.6° 106° SIE» 7/57—12/58 
Tees tes. sees +53. 49° nO Bey Wie | Sees? 83° 100 7/57—12/58 
WEWISONY, 650006 Guide —67° 36° 622586. ESS 104° SE 7/57—12/58 
Murchison Bay .... + 80° 03’ TOME 0 +72.2° TI Sule 9/57—12/58 
Oftawarn ern. +45° 24’ 75° 54 W | +56.8° SELS 100 7/57—12/58 
SUPREME sdseac dE 127 115° 36’ W | +58.2° 300° 2283 7/57—12/58 
Berkeley. ci. ce 5705200) 1222 Sew 4470 298° Sale 7/57— 9/58 
IResolute nr re +74° 41’ 94° 54’ W | +82.9° 289° SL. 7157 -12/58 


to nine stations for which I.G.Y. records are 
relatively free of breaks in continuity and for 
which the counting rates of the neutron moni- 
tors ensure reasonably small statistical errors. 
Table I lists the stations selected, together 
with their geographic and geomagnetic coor- 
dinates. 

By regarding as permissible the interpolation 
of no more than two missing bihourly figures 
to complete at any station an otherwise 
acceptable day, the total number of days of 
simultaneously available records found in the 
eighteen separate months ranges from 18 to 31, 
yielding satisfactory statistical accuracy. Inter- 
polations have been made where necessary by 
fitting values to the general trend shown by 
adjacent bihourly figures. 


Harmonic analyses carried out 


The analyses for each station have been 
carried out using the facilities of the BESK 
electronic computer in Stockholm. The basic 
data fed to the computer for each analysis 
consisted of thirteen consecutive bihourly count 
totals. The thirteenth value was obtained by 
considering each day to extend from oo hrs 
to 24 + 2 hours. The computer automatically 
corrected for the difference in the ıst and 
13th values on the assumption that it represents 
a linear secular intensity change over the 
mean day. The programme was arranged to 
yield the amplitude in per cent and the times 
of maximum of both the diurnal and semi- 
diurnal components, the usual procedure of 
least squares fitting being followed. All anal- 
yses were based on days commencing at 00 
hrs Greenwich Mean Time, and used pressure- 
corrected data. 


Two main sets of analyses were made; the 
first used data from all available days of simul- 
taneous records and will be referred to as 
the “all selected days” analysis, while the 
second used data from only those selected days 
remaining after exclusion of “disturbed” days 
and will be called the “quiet selected days” 
analysis. The classification of days as disturbed 
has been based on visual inspection of plotted 
records and not on any strict criterion. In 
general the disturbed days are those including 
the onset and decreasing phase of Forbush- 
type decreases and perhaps one or two following 
days, together with any other days on which 
marked anomalous primary intensity changes 
occurred. 

The grand totals for each station over the 
whole period have also been analysed both 
for “all selected days” and “quiet selected 
days” as have also the totals for the 46 days 
classified as disturbed. 

Figures I and 2 present the results for the 
diurnal component plotted in the conventional 
manner on vector summation dials for the 
“all selected days” and “quiet selected days” 
cases respectively. For each station the whole 
vector summation plot has been rotated clock- 
wise through an angle which brings all plots 
approximately parallel and in the general di- 
rection of 12 hrs G.M.T. on the dial. This 
facilitates comparison of phase changes at the 
different stations. These angles of rotation are 
as follows. 


Hobart 99° Murchison Bay 335° 
Resolute 175° Mawson 338° 
Berkeley 182° Leeds 344° 
Sulphur Mt. 188° Uppsala 346° 
Ottawa ee ie 
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MURCHISON 
RESOLUTE BAY 


SULPHUR MT 
MAWSON 


HOBART 


All vectors for each station have been 
rotated through fixed angles (see text) 


Hig. T- 


453 


UPPSALA 
OTTAWA 


BERKELEY 
LEEDS 


Vector summation diagrams for the observed monthly mean diurnal 


variation in the nucleonic component intensity at 9 stations during the I.G.Y.: 
“all selected days’’ analysis. The month numbers identify corresponding vectors 
at different stations. 


3. The expected phase differences between 
stations. 


It is evident in both Figures 1 and 2 that 
the month to month phase changes at the 
various stations are not well correlated. The 
changes are less irregular in Figure 2 than in 
Figure 1 indicating that the exclusion of 
disturbed days improves the general agree- 
ment, even though purely statistical errors are 
made relatively larger, but the detailed changes 
are still not very well correlated. However it 
is clear that over the whole 18 month period 
no really pronounced semipermanent changes 
in direction of the primary anisotropy occurred. 

In order to judge the measure of phase 
agreement between stations in any one month 
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we must establish the “expected” or “standard” 
phase differences in Universal Time between 
the stations. If the daily variation were due 
entirely to an anisotropy in the primary radia- 
tion entering the geomagnetic field, then 
apart from the uncertainty introduced by 
statistical errors, we would expect to observe 
quite specific fixed intervals between the 
intensity maxima at different stations. In 
principle we could calculate these intervals 
from a knowledge of the energy dependence 
and spatial form of the anisotropy, the trajec- 
tories of primaries traversing the earth’s field, 
the complete directional sensitivity pattern at 
the top of the atmosphere for each recorder, 
and the geographic coordinates of the recording 
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MURCHISON 
RESOLUTE BAY 


SULPHUR MT. 
MAWSON 


HOBART 


All vectors for each station have been 


rotated through fixed angles (see text) 
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LEEDS 
BERKELEY 


OTTAWA 
UPPSALA 


Fig. 2. Vector summation diagrams for the observed monthly mean 
diurnal variation: “quiet selected days” analysis. 


stations. However as we can only infer the 
properties of the anisotropy from our ground 
level observations, it seems more satisfactory 
and practical to assume for our present purpose 
that the “expected” U.T. intervals are those 
exhibited by the mean daily variations over a 
relatively long period. 

Table 2 A sets out the phase angles of the 
long term (12 and 18 month) mean diurnal 
variations determined by different methods 
(see description with the Table). The phase 
angles quoted are values of fi appearing in 
the expression R, sin (O+f,) for the diurnal 
component, with the time zero, © =o, at oo hrs 
G.M.T. Table 2B lists the corresponding 
differences in phase angle between stations 
arranged approximately in the order in which 


the maxima occur (15°=1hour). These differ- 
ences are quite similar for the various cases. 
However in order to minimize any effects of 
possible annual or seasonal variation in ampli- 
tude and the effects of inclusion of disturbed 
days, we have chosen to determine the 
“standard” phase differences from the two 
overlapping 12-month quiet day periods 
(columns g and h) for which we have strictly 
comparable records for all 9 stations. The phase 
determinations in these cases have been made 
directly from the vector plots of Figure 2, as 
the phase angle of the resultant, or mean, of 
the appropriate 12 vectors. An average is 
chosen from the two 12-month periods and 
the accepted standard phase differences are 
given in the final column¥of Table 2 B. _ 
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Table 2. 


2A Phase angles of long term mean diurnal variations in the nucleonic component. 
(Phase angle /,°, in the expression R, sin (0 + /,) for the diurnal component with 0 = o at oo G.M.T. 
2 B Differences in phase angle. These are reckoned positive clockwise on a U.T. dial. 


a) From analysis of grand totals for 18 months , all selected days 
b) » » » » » » 18 » , quiet selected days 
c) From mean of 18 monthly vectors , all selected days 
d) » » » 18 » » 3 , quiet selected days 
e) » > » 12 » » , 9/57—8/58, all selected days 

f) » ER 2 » » » 10/57—9/58,° » » » 

g) » by he ans » » » 9/57—8/58, quiet selected days 
h) » » » I2 » » » 10/57—9/58, » » » 


i) The adopted standard phase differences (approximate means from columns g and h) 


2A | a | b Cc | d e | f | g N h | i 
Hobaressen ger... 22 15 23 16 19 16 10 8 
ippsalas.c. << ares « 256 256 257 256 259 258 257 255 
Léeds shoes o's sce ices s 255 256 256 255 255 259 253 254 
MAWSON. essen 245 245 247 247 252 249 248 248 
Murchison Bay 251 243 250 240 
Ottawa ER asie 154 153 154 153 153 149 153 149 
Sulphur Mt... =< 104 IOI 106 103 106 97 102 94 
Berkeley re... 97 95 93 OI 
PRESONICG 6 core. re seu 96 86 99 89 98 87 88 83 
25 
EI URN rn. 2,02 126 119 126 120 120 118 113 113 113 
CEE ae eee I o I I 4 —I 4 I 2 
ae Mer re. Es. Le 10 ET 9 8 3 Io 5 6 6 
MEME eds. „on | Sr N en } 55 } 3 I 6 —2 8 3 
BEE ee u à 98 94 97 91 94 
Re ecke 50 52 48 50 47 52 51 55 53 
Benannte 9 2 9 3 6 
ky en... } : } 15 } 7 } + 23 8 5 8 7 
HR Malen 74 71 76 73 79 71 78 75 76 
sm BR Figure 3 shows the standard phase intervals 
on a 24-hour U.T. dial. We are at present 
interested only in the phase differences, i.e. 
the relative and not the absolute positions of 
. the vectors, and thus the time scale is arbitrary. 
Wi ees Sasa H 4. Determination of a U.T. component of 
) the diurnal variation. 
Having established expected phase differ- 
H Hobart ences in U.T. between the stations, we can 
x N proceed as mentioned in Section I to see 
M Mawson whether for each month separately, we can 
MB Murchison Bay improve the agreement between the observed 
Qs0teve and expected phase intervals. 
P P 
MB SM Sulphur Mt. 1 \ J 
MLU B Berkeley We commence with the hypothesis that 
R Resolute 


Fig. 3. 24-hour dial showing the adopted “standard” or 

expected phase relationships between the diurnal varia- 

tions at different stations in Universal Time. The phase 

differences in angular measure are given in the final 
column of Table 2. 
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the monthly mean diurnal variation contains 
a component arising from synchronized varia- 
tions in U.T. affecting all stations. The maxi- 
mum. of this component would thus occur 
at the same U.T. at all stations. We make the 
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simplifying assumption that it has the same 
amplitude for all stations considered. This 
seems justified in view of the fact that such 
transient and approximately simultaneous 
world-wide primary intensity changes as For- 
bush-type decreases produce almost the same 
percentage changes at all stations above the 
knee of the latitude curve. The possibility of 
making this simple assumption has led to our 
choice only of relatively high latitude stations 
for the investigation. 

We shall examine the results of the “all 
selected days” analysis. For a particular month 
the nine station diurnal vectors are drawn 
on a single U.T. dial on transparent paper. 
This is then placed over the standard dial of 
Figure 3 so that the origin points coincide. 
In most months quite marked differences from 
the standard phase intervals are evident and 
the fit is poor no matter into what relative 
positions the two dials are rotated. However 
by such relative rotation about the origin 
point it is possible in general to find a position 
where subtraction of a single fixed vector of 
suitable amplitude from each of the experi- 
mental ones will alter their phase angles in a 
manner resulting in considerably improved 
agreement with the standard intervals. Although 
the method is rather crude, the appropriate 
vector and orientation of the dial giving best 
improvement are not difficult to find, and 
more elaborate fitting methods are not war- 
ranted. We call the vector found the U.T. 
vector for the month considered. (It applies 
of course only to the mean diurnal variation 
for the limited number of days used in this 
analysis and cannot be considered applicable 
for any other selection of days from the 
month.) It should be pointed out that if we 


Table 3. The monthly mean U.T. diurnal comp- 
onents determined for the ‘‘all selected days” 
analysis, R,sin(0+-f,) with 0 — o at 00G.M.T. 


Month | nl tae | Month | RO || Vie 
7/57 .100 120 4/58 | ..050 300 
8 .100 180 5 2120 fe) 
9 .200 180 6 oO — 
Io .100 140 7 .080 32 
II 050 90 8 .250 125 
12 fe) — 9 -100 65 
1/58 150 305 10 .200 125 
2 025 150 Il .080 195 
3 100 75 12 .100 100 
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Fig. 4. Summation diagram for the monthly U.T. 
vectors the subtraction of which produces the most 
improved inter-station phase agreement. 


were considering only two or three stations 
it would be possible to determine a unique 
vector which would produce exact agreement. 
For more than three stations we can find a 
vector which on subtraction gives only the 
most improved general agreement. 

Details of the U.T. vectors found are given 
in Table 3 and the vectors are plotted on a 
summation dial in Figure 4. Although the 
accuracy of the method used is not great, it is 
clear that in many cases the amplitude of the 
U.T. vector is relatively large and its sub- 
traction can therefore make very appreciable 
changes to both the relative amplitudes and 
phases of the original vectors. An example 
illustrating these effects for the month August 
1958 is shown in Fig. 5. Figure 5A shows 
the original vectors together with the U.T. 
vector to be subtracted. This is a case in 
which the original phase intervals are greatly 
different from those of the standard dial in 
Fig. 3. Fig. 5 B shows the resultant vectors, 
and the direction marked H, is that for the 
Hobart vector on the standard dial. It indicates 
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GMT 
00 


SM 


Fig. 5. Harmonic dials for August 1958, illustrating the 
possible large effects of subtraction of a U.T. component 
common to all stations. The station vectors are identified 
as in Fig. 3. 5 A: original vectors. 5 B: after subtraction 
of the vector V. The direction marked Hs (=Hobart 
standard) indicates the orientation of the standard dial 
of Fig. 3, giving best phase fit to the resultant vectors. 


MURCHISON SULPHUR MT. 
RESOLUTE BAY = Mawson HOBART 
1 ; 3 
7 
6 
6 
12 
GMT. 
00 
18 06 
0.2% 


All vectors for each station have been 
rototed through fixed angles (see text) 


OTTAWA 
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the orientation of the standard dial which 
gives best fit to the corrected vectors. The 
fit in this case is still far from perfect but a 
considerable improvement has certainly been 
made. 

Fig. 6 shows the summation dials for the 
“all selected days” analysis after correction for 
the effects of the U.T. vectors, and is to be 
compared with Fig. 1. It is clear that the 
month to month phase changes are now 
much more consistent at the various stations 
as one would expect if the vectors now 
describe much more accurately the diurnal 
component arising from a primary anisotropy 
alone. 

The degree of consistency in Fig. 6 is even 
more satisfactory than that in Fig. 2 (for the 
“quiet selected days”). By following a similar 
procedure month by month for the quiet day 
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Fig. 6. Vector summation diagrams for the monthly mean diurnal variation: 

for the “all selected days’’ analysis, after subtraction of monthly Universal Time 

components common to all stations. (Note that in July 1957 the Resolute vector 
has almost zero amplitude.) 
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case, U.T. vectors of appreciable amplitude, 
often similar to those for the “all selected 
days” case can still be found but the statistical 
accuracy is somewhat less. Thus the U.T. 
component is not a feature of the disturbed 


days only. 


5. Changes in direction of the primary 
anisotropy. ; 


If the U.T. component found is a genuine 
effect which shifts the phase of the diurnal 
variation arising from a primary anisotropy 
alone, then the positions of the corrected 
vectors on the U.T. dial and not those of 
the original vectors will be directly related to 
the direction of the primary anisotropy. 
Further, the month to month changes in the 
optimum orientation of the standard vector 
dial with respect to U.T. will determine the 
changes in direction of anisotropy with which 
the 9 stations best agree. That is, we obtain an 
average determination based on records from 
9 stations. 

We may take for example changes in the 
phase angle of the standard Hobart vector as 
our best estimate of changes in the direction 
of anisotropy. The month by month phase 
behaviour of this vector is shown in Fig. 7, 
which has been drawn as a summation dial 
with all vectors given an arbitrary fixed 
amplitude. Except for a constant additive 
angle the plot of Fig. 7 shows directly the 
changes in the effective direction of anisotropy. 
Clockwise changes in vector direction indicate 
that the direction of maximum intensity has 
shifted further from the earth-sun line in the 
sense of the earth’s rotation. 


Fig. 7. Showing the month by month changes in direction 
of the primary anisotropy, deduced as described in Section 
5. The actual direction of the anisotropy in space relative 
to the earth-sun line would be obtained by adding a 
constant angle to the phase angle of each vector. A clock- 
wise change in vector direction corresponds to a shift 
in the direction of anisotropy in the sense of the earth’s 
rotation. 


We conclude from this Figure that the 
direction of maximum. intensity remained 
virtually constant from July to December 
1957, then changed by about 25° for the 
next few months before becoming rather 
variable until the end of 1958, although the 
mean direction for the latter half of 1958 was 
approximately the same as that for the similar 
period in 1957. Attention should be drawn 
here to Figure>4 which shows that marked 
changes in the phase trend of the U.T. vectors 
also occurred at about the same times as 
changes in the direction of anisotropy (Fig. 7). 
It is also interesting to notice that December 
1957 coincided approximately with the lowest 
primary intensity reached during the present 
solar cycle. The existence of some degree of 
quasi-persistent behaviour in the U.T. vectors 
as shown in Figure 4 and the possible associations 
with the direction of anisotropy and the mean 
level of primary intensity may be significant, 
but the present evidence does not justify further 
speculation here. 

Although we have not been concerned with 
the absolute direction of the primary anisotropy 
nor with the actual mean deflections in geo- 
graphic longitude suffered by the primaries 
affecting each station, it is nevertheless of 
interest to determine the relative deflection 
angles. For instance if we represent the deflec- 
tion angle measured eastwards for Hobart 
by Y, then from our standard phase intervals 
and the geographic longitude of the stations, 
we deduce that the effective mean deflection 
angles are as follows. 


Resolttés «eed tink ee de Y —42° 
Mawsairwu.. di er cae Y — 39% 
Sulphur Me fe OPA RUES PS 
Beckeloys Al aa ame tes oe Y- 8° 
Hobart. enic? a th er of 

Ottawa cl. hu Re ETS W+ 5° 
Murchison Bayer Gye: ee er Pr 
Uppsala steh een a ns D +16° 
Leeds ae dust lat D +34° 


These estimates are of course based on the 
observed mean time of maximum of onl 
the diurnal component and will remain valid 
only if the energy dependence of the primary 
anisotropy remains substantially constant. It 
will be interesting to compare these relative 
values with theoretical predictions when these 
become available. 
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6. Origin of the U.T. component. 


The apparent existence of a diurnal variation 
synchronized in U.T. at all stations does not 
necessarily imply that the primary radiation 
undergoes world-wide pulsations with a true 
solar-day periodicity. There are few ways in 
which world-wide modulation of primary 
intensity could conceivably be controlled 
directly by the earth’s rotation—one could 
perhaps envisage that the rotation of the 
eccentric geomagnetic dipole field in a nonuni- 
form medium might produce small periodic 
changes in a geocentric electric field. However 
it is unlikely that any such effects would 
be consistent with the large phase changes 
shown by the diurnal U.T. variations found 
(Table 3). 


It must be remembered that in determining 
the expected phase intervals due to a primary 
anisotropy alone, we made the implicit assump- 
tion that any disturbing effects would approxi- 
mately cancel in an annual mean for quiet 
days. If there were any tendency over a long 
period for a U.T. component to exhibit a 
preferred phase, this assumption would be 
unjustified. It would therefore be valuable to 
have theoretical estimates of the phase inter- 
vals to compare with the experimental values 
used here. Mc Cracken (private communica- 
tion) has made estimates for some stations 
which do in fact agree fairly well with the 
intervals used here. For instance he has pre- 
dicted (Mc CRACKEN, 1959) that although 
Mawson is 45° (3 hours) east of Uppsala, the 
latter should record a maximum in the diurnal 
variation of the nucleonic component some 20 
minutes (5°) earlier in U.T. than Mawson. 
This compares favourably with the lead of 8° 
used in the present investigation. 


Thus it seems likely that the U.T. component 
exhibits no long term persistent features, and in 
a monthly period may arise from irregular 
primary intensity changes unrelated to the 
earth’s rotation. We shall briefly examine some 
possibilities. 

The pressure corrected bihourly intensity I (t) 
measured by a neutron monitor during a 
day may be represented as the sum of several 
time-variable functions, e.g., 


I(t) =I+a(t)+b() +e() +d (0) 
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a(t) is a function containing all harmonics 
describing the daily variation arising 
purely from a primary anisotropy. It is 
essentially a local time phenomenon. 

is a function describing purely statistical 
fluctuations, independent at different 
stations. 

c(t) describes the contribution made by rela- 
tively long term secular changes (e.g. 
the 27 day and month to month changes 
in mean intensity level). 

describes non-random but irregular short 
term variations of a quasi-persistent char- 
acter, leading to serial correlation of 
consecutive bihourly intensity values. 
These may be due to rather localized 
effects or to synchronized world-wide 
variations (e.g. small Forbush decreases and 


solar flare effects would be of this type). 
Only the functions c(t) and d(t) need be 


considered as possible sources of an apparent 

U.T. diurnal variation in a monthly average. 

The function c(t) could produce such a com- 

ponent as a result of the general mean curvature 

of changes of intensity level occurring in a 27- 

day type of variation. This type of “curvature 

fect” has been discussed by CHAPMAN and 

BARTELS (1951). Depending on its sense on the 

majority of days in the period considered, 

the general curvature will introduce a com-. 
ponent to the mean daily variation with 

either maximum or minimum at the centre 

of the daily period chosen. Moreover the 

phase of this component will be independent 

of the choice of starting point for the daily 

intervals. Since the 27-day variations are 

world-wide this effect would produce a diurnal 

variation synchronized in U.T. at all stations. 

The magnitude of this curvature effect may be 

determined simply in any particular case. In- 

spection of a plot of daily mean intensity 
values shows that the month in our analysis 

likely to contain the greatest curvature effect is 

August 1957, although it is far from an ideal 

case. In this month actual calculation of the 

amplitude of the diurnal component intro- 
duced by the curvature effect shows it to be 

less than 0.005 per cent. 

Since the curvature effect is negligible, we are 
left then with only the function d(t) as a 
possible source of our U.T. component. That 
approximately synchronized world-wide in- 
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tensity variations of this type occur at least 
on some days is certain. Forbush decreases are 
obvious examples, but other less prominent 
cases may be readily seen when data plots for 
many stations are compared. Quite small varia- 
tions of a world-wide character are difficult 
to recognize since they would be partly 
masked by the local time dependent and purely 
random fluctuations present simultaneously. 
Nevertheless it seems reasonable to assume 
that a high proportion of days will be affected 
to some extent in this way. This assumption 
is strongly supported by the fact that signif- 
icant and highly correlated changes occur at 
different stations in daily mean intensity values. 
In order to estimate the possible amplitude of 
the U.T. diurnal component arising from 
this type of intensity fluctuation we may 
consider an artificial model which nevertheless 
simulates quite well the typical pattern of 
quasi-persistent variations actually observed. 


Suppose we have a long time series con- 
sisting of equally spaced independent values 
of a random normal variable whose standard 
deviation is o, and that we then form from 
it a new series consisting of the running sums 
of 4 consecutive values. The new series will 
have standard deviation 20 but will possess 
the property of quasi-persistence, with each 
value being independent only of those four 
or more on either side of it. If o is suitably 
chosen, trials show that such a series does in 
fact resemble quite strikingly the type of 
irregular fluctuations actually observed in 
bihourly cosmic ray intensity figures, although 
in the cosmic ray case the range of the fluctua- 
tions may change considerably from day to day. 


If we take independent samples of 12 con- 
secutive values from the original random series 
and fit a first harmonic to each sample, it can be 
shown that its root mean square amplitude 
will be 0.577 o (Ch. 16, CHAPMAN and Bar TELs, 
1951). If we now repeat this procedure using 
similar samples from the derived quasi-per- 
sistent series the root mean square amplitude 
will be 1.752 o. Suppose now that we exam- 
ine the effect in the mean over N samples 
(corresponding to N days in the cosmic ray 
case). In the case of the original random series 
the root mean square amplitude in repeated 


selection of N samples will then be 0.577 o/Y N. 
In the case of the quasi-persistent series the 
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expected amplitude will depend on whether 
the N samples are consecutive in the series 
(in which case they are not independent), or 
whether they are independent and not adjacent 
in the series. (In the cosmic ray case there are 
usually some missing days in an otherwise 
consecutive group). If the N samples are 
independent the root mean square amplitude 
in repeated selection of groups of N will 


be 1.752 o/VN, whereas if all the samples are 
consecutive the value‘is slightly higher, and is 
given by 


1.414 (1.866 N - 0.331)? o/N. 


To relate these figures to the cosmic ray 
case we must choose a value for o which 
yields quasi-persistent fluctuations of a rea- 
sonable magnitude not inconsistent with the 
observed type of intensity variations. 0—=0.3 % 
is such a reasonable choice. It is in fact less 
than the statistical standard error of bihourly 
intensity values (~0.4 %) for typical record- 
ers considered in this paper. For a typical 
“monthly” mean we may take N=25. These 
values would yield a root mean square ampli- 
tude for a fitted diurnal variation of between 
0.105 % and 0.115 % (for the non-consecutive 
and consecutive-day cases respectively). 

Obviously these values would become still 
larger if either o or the persistence period were 
increased slightly. However the example 
chosen is sufficient to show clearly the rela- 
tively large contribution which such variations 
could make to the observed daily intensity 
variation. Moreover since we are considering 
them to represent approximately synchronized 
intensity variations at all stations their contri- 
bution to the diurnal variation would have 
the same time of maximum in U.T. at all 
stations, and the phase would naturally be 
subject to large erratic changes from month 
to month. As the U.T. components found in 
our analysis (Section 4) have an average ampli- 
tude of about 0.1 % it seems that they could 
be explained quite well by irregular short-term 
modulation of the world-wide intensity. 


7. Conclusion. 


The evidence presented suggests that a 
Universal Time component plays an impor- 
tant role in determining the form of the 
observed daily intensity variation. Subtraction 
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of such a component from monthly mean 
diurnal variations observed at different stations 
leads to more consistent month to month 
phase behaviour at the various stations and 
considerably improved interstation agreement 
on the direction of the primary anisotropy 
responsible for the principle component of 
the daily variation. 

It has been shown that the U.T. component 
may result from the existence of irregular 
world-wide primary intensity changes of the 
short term quasi-persistent type and it may 
not be necessary to postulate daily periodic 
modulation of the primary intensity caused 
directly by the earth’s rotation. Nevertheless 
examination of further experimental data in 
the light of these suggestions is required. 

We have considered only the diurnal varia- 
tions. In principle we could follow an exactly 
similar procedure for the semidiurnal or 
higher harmonics since each harmonic is 
formally independent of the others. However 
statistical uncertainties are too great in these 
cases to make this worth while. It should 
however be pointed out that the proposed 
origin of the U.T. variation would produce a 
semidiurnal component of mean amplitude 
similar to that for the diurnal component. 
Its effects would therefore be relatively greater 
on the phase and amplitude of the smaller 
semidiurnal variation (see Fig. 8) and may 
provide an explanation of the very erratic 
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Fig. 8. Vector summation diagrams for the observed 
monthly mean semidiurnal variation in nucleonic com- 
ponent intensity: “all selected days’’ analysis. Consecutive 
vectors are for months 7/57—12/58 for all stations except 
Berkeley (7/57—9/58), and Murchison Bay (9/57—12/58) 
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and inconsistent behaviour of this variation at 
different stations. 

In a similar way the presence of short term 
irregular world-wide intensity variations un- 
related to any primary anisotropy could 
provide at least a partial explanation of the 
high degree of day to day variability in the 
observed daily variation and the difference in 
behaviour at different stations. Such an ex- 
planation could remove some of the difficul- 
ties inherent in attempts to explain the variabil- 
ity in terms of rapid changes in an anisotropy 
of complex form. 

Other features of the daily variation in 
which inconsistent behaviour has been ob- 
served are the changes related to the degree 
of geomagnetic activity. In the majority of 
published results, increased amplitudes and 
earlier times of maximum are reported to 
accompany increased geomagnetic disturbance 
but the reverse is occasionally observed at 
some stations. In the present analysis for in- 
stance, on comparing the 18-month mean 
diurnal variations for “quiet” days with the 
mean variations on the 46 days classified as 
disturbed (excluding Murchison Bay and 
Berkeley for which we have less than 18 
months data) we find the differences to be 
rather inconsistent at the various stations. 
Although the disturbed days were not selected 
on a magnetic disturbance criterion they 
nevertheless do correspond in almost all 
cases with enhanced geomagnetic activity. By 
following the same procedure with the mean 
diurnal components for the disturbed days as 
that used for monthly means, a U.T. compo- 
nent of amplitude 0.20 per cent and suitable 
phase has been found which on subtraction 
from all station vectors produces a very 
marked improvement in agreement with the 
standard phase intervals between stations. 
Whereas compared with the corresponding 
quiet day variations the original disturbed 
day diurnal components show maxima ranging 
from 3 hours earlier (Hobart) to 0.5 hours 
later (Leeds) with only small insignificant 
differences for Uppsala, Mawson, and Ottawa, 
the resultant vectors after subtraction of the 
U.T. component all exhibit maxima between 
1.1 and 2.1 hours earlier. 

The various effects discussed show clearly 
that the possible influences of a U.T. compo- 
nent chould not be overlooked in any interpret- 
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ation of experimental observations of the daily 
intensity variation. In particular, care should 
be exercised when conclusions about the direc- 
tion of primary anisotropy are drawn from 
the phase of the daily variation observed at 
ane stations. 

It must be stressed that the reality, relative 
importance and origin of the component have 
not been firmly established, and although the 
success of the hypothesis is strongly suggestive, 
independent methods of testing it are required. 
From the discussion of Section 6 it is clear 
that if world-wide short-term intensity fluc- 
tuations occur, even if small and somewhat 
intermittent, they must give rise to a com- 
ponent of the daily variation. They undoubt- 
edly occur on some days, but could well 
occur on many others, partly masked by 
local time dependent and statistical fluctua- 
tions. Perhaps the best method of separating 
them for study would be to average intensity 
data in Universal Time for a number of 
neutron monitor stations suitably situated at 
roughly the same latitude but different longi- 
tudes so that the averaging process would 
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eliminate the local-time variations caused by a 
primary anisotropy. 
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NOTE 


The detailed numerical results of the harmonic anal- 
yses described in Section 2, together with the statistical 
standard errors of estimate of amplitudes and times of 
maximum can be made available in tabulated form to 
those wishing to make use of them. 
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ERTTERS- TO VHEPEDITOR 


“The Yearly Circulation of Chloride and Sulfur in 
Naturel, Meteorological, Geochemical and Pedological 
Implications”, by Erik Eriksson, Tellus, 11, PP. 
375—403, and 12, pp. 63—109 (1959—1960). 


Dear Sir, 

In his paper Eriksson makes a worthwhile attempt 
to combine a large number of observations about 
sea salt components into a comprehensive picture. 
Included in his discussion is some of our data from 
the U.S. rain water analyses project. In a number of 
cases, the conclusions drawn in this paper differ 
from those of the original authors. Most of this 
disagreement originates because Eriksson prefers 
to present these data in a different way. Since this 
seems to touch on a more fundamental question in 
the interpretation of rainfall data, a few remarks 
may be justified. 

Our interpretation of rain analysis data (JUNGE 
and WERBY 1958) was based on maps of average 
concentration of various ions, (Cl-, Nat, K+, SO, — 
etc.), in rain, expressed in milligrams per liter. The 
purpose was to obtain from such data, if possible, 
information about the origin and source areas of 
these atmospheric constituents. Eriksson uses for 
the same purpose maps of total deposition of these 
constituents, expressed in kilograms per hectar 
per year. These latter maps are calculated from the 
concentration maps by multiplication with the 
average yearly amount of rainfall. As the title of 
his paper indicates, he is concerned with land-ocean 
circulation of such constituents, and, for yearly 
budget considerations, maps of the total annual 
deposit are undoubtedly of great value. If they are 
used, however, for interpretation of the origin of 
these components, objections must be raised. 

The main reason that we used rain water concen- 
tration maps for our interpretation, is that they 
reflect the average tropospheric air concentration 
much better than maps of total deposit. Although 
we gave specific reasons for this in the publications 
of our rain water results (JUNGE and WERBY 1958, 
JUNGE and GUSTAFSON 1957, JUNGE 1958), Eriksson 
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restricts his discussion of this question to the sen- 
tence: “The concentration is not a good indicator 
because it is strongly dependent upon the amount 
of water vapor present in the air, subject to conden- 
sation” (page 56). We would like to discuss this 
point in more detail. The amount of water available 
for the removal of the chemical constituents from 
a certain volume of cloud air is largely controlled 
by the liquid water content of this cloud at the 
onset of rain. This value varies with the type of 
cloud, but in the rather narrow limit of about 1.5 
to 3.0 grams per cubic meter. In warmer climates, 
where the water vapor content of the air is higher, 
this liquid water content is released by a smaller 
amount of cooling than in colder areas. There is a 
small, but systematic decrease in the average liquid 
water content of precipitating clouds, from the 
equator to the poles, because of a variation in the 
composition of the cloud population. In warmer 
climates, there is generally a higher percentage of 
convective clouds which tend to have a higher 
liquid water content. But this is a smooth and 
gradual decrease, with a gradient much smaller than 
those exhibited on most of the concentration maps 
of rain water constituents. 

Basically, the material which is contained in rain 
water comes from two predominant sources: 
Collection of material within the cloud, during 
the formation process of the cloud and rain drops, 
and secondly, from the interception of aerosol 
particles by falling rain drops below the cloud. Unless 
the original air concentration of the aerosol is much 
larger below the cloud than within the cloud (i.e., 
sea salt particles over the ocean and along the coasts), 
the first source is the predominant one (JUNGE and 
GUSTAFSON 1957). This can be expressed quantita- 
tively. If, within the cloud, / is the liquid water 
content (gr/m?), C the concentration of the chemical 
constituent in air (g/m), Kthe same in rain (g/cm? 
=mg/liter) and & the fraction of material removed 
from the air within the cloud, we obtain: 


ar 
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During the periods of months or longer, the fluc- 
tuations of &/l with cloud type, etc., will average 
out, and (3) will vary only slightly with location, 
season and surface (land, ocean). As a good approxi- 
mation, therefore, K will be proportional to C. If, 
for sea spray, C =10 ug/m*, ex 1.0 and /=2 g/m?,we 
we obtain K=s mg/liter, a figure which is of the 
right order. The increase of this concentration by 
interception of particles below the cloud—if im- 
portant—will, again, on the average, vary only 
slightly with latitude, season and surface (land, 
ocean). 


The question now arises, to what extent the air 
concentration is modified by preceding rainfalls 
and washout or by the average amount of rainfall 
if we think in terms of geographical distribution 
rather than air mass history. There is no doubt 
that in an individual precipitating cloud or cloud 
system, constituents can be almost completely re- 
moved. This is demonstrated by the fact that during 
individual rains, the concentration decreases con- 
siderably with time. But if we deal with large air 
masses, as we do when we use monthly averages, 
the distribution of rain in time and space is such that 
only a small fraction of the constituent is removed per 
day. With ¢ =1 we calculated residence times (JUNGE 
and GUSTAFSON, 1957) of 6 days for an average 
rainfall intensity of 0.1 inch per day, and a liquid 
water content of 2.0 g/m, conditions which pre- 
vail over most of the U.S. These calculations with a 
simplified model give the right magnitude. Obser- 
vations indicate a range from $ to 40 days, which 
can be explained by an & smaller than 1, and by 
simplifying the assumption for other parameters. 
Residence times of this order show that, on the 
average, even when air masses cross large areas like 
the U.S., the constituent in the air is not decreased 
considerably. 


Let us now consider two extreme cases of washout 
in the atmosphere, and the effect on concentration 
maps as well as on total deposition maps. We assume 
a certain constant geographical distribution of con- 
tinuous sources of chemical constituents at the 
ground. In the first case, the residence time of the 
pollutant may be extremely short, so that it will 
be removed before it can spread out to any degree. 
It is obvious that in this case the concentration, as 
well as the total deposition, will reflect very closely 
the distribution of the sources. Both maps could 
be used for source interpretation. In the second 
case, the residence time shall be long, so that con- 
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siderable horizontal mixing and transport occurs 
before any appreciable amounts are removed. The 
air and rain concentrations will now show almost 
uniform values, and the total deposit will reflect 
the map of total rainfall. In reality, conditions are 
somewhere between these extremes. With residence 
times of 6 to 40 days, it should be expected that 
the air concentration gives only a blurred picture of 
the source distribution, with extension of the areas 
of main concentfation downwind. But the maxi- 
mum concentration will still be close to the source 
area, and the pattern will give information about 
the main source and sink areas. For total deposit, 
however, the pattern of the rainfall distribution is 
superimposed upon that of air concentration, and a 
reliable interpretation becomes more or less im- 
possible. 

The example of Sr9° may illustrate this. It is 
now commonly agreed that, over wide areas of 
uniform climatological structure, the deposition of 
Sr°0 is proportional to the total amount of rainfall. 
The average release of Sr®° into the troposphere is 
primarily latitude dependent, and because of the 
speed of circulation and of residence times of 4 
weeks, the air concentration is fairly uniform within 
latitude belts. The higher rainfalls along the west 
coast of the U.S. result, therefore, in higher Sr°° 
deposition than further inland, but nobody would 
draw the conclusion that Sr?® comes from the sea. 

In a similar way the conclusions Eriksson draws 
from total rainfall deposition maps as to the origin 
of the constituents are doubtful in some cases, and 
incorrect in others. 

His map 7.5 of the excess sulfur in precipitation 
over the U.S., exaggerates the influence of indus- 
trial sulfur and obscures the gradients along the 
coast. Because total rainfall changes considerably 
with the distance from the coast, we cannot expect 
that this map can bear out our conclusion that most 
of the sulfur is of continental origin. The concen- 
trations of excess sulfur, however, are clearly lower 
over the oceans than over the continents. 

Similarly, we cannot agree with the conclusion 
(pages 62 and 64) that excess Nat and K+ come 
from the sea (pages 62 and 64), because we think 
that his maps 8.2 and 8.4 (to stay just with the U.S. 
data) should not be used for this argumentation. 
In addition, our data on which these maps are based, 
are not regarded reliable enough along the coasts 
to draw any conclusions from them. We indicated 
(Junge and Wersy 1958, page 420), that our data 
for the excess Nat and K+ along the coasts are not 
too reliable because they appear as small differences 
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of two comparatively large figures, whereas they 
become very reliable further inland. A similar 
caution should be exercised with most of the other 
excess data over or near the ocean. For this reason, 
we did not use the pattern of the excess concentration 
maps of Nat and K+ for our considerations. Our 
line of reasoning was based on the important fact 
that the average concentrations of Nat and Kt 
over most of the U.S. are very much the same despite 
the fact that the ratio K+/Nat is only 0.041 in sea 
water and thus 20 times smaller than in rain. There 
is no basis for the assumption that the ionic composi- 
tion in sea spray particles deviates that much from 
sea water. On the contrary, coastal rains in unpolluted 
areas have ratios close to that of sea water. We 
concluded, therefore, that most of the K+ must 
come from the soil. Because of the similarity in 
chemical behavior of Nat and K+, we suggested 
that it is very likely that the excess Nat is also 
derived from the soil. 

As we pointed out earlier, KALLE (1953—1954) 
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in a very interesting study, provided considerable 
evidence of the mineral origin of the excess Nat. 
K+, Cat+ and other cations in rain over Europe, 
It would have been very useful if this paper had been 
included in Eriksson’s discussion. Also, the reader 
would have appreciated a reference of METHAM’S 
(1950) paper on the SO: concentration and SO, 
deposition over England in connection with Figures 
7-7 to 7.10. Metham’s work indicates that industry 
is an important, if not the predominant source of 
sulfur over England, in contrast to the conclusions 
drawn from Figures 7.7. to 7.10. The heavy rain- 
falls on the west coast of Ireland certainly increase 
the deposit, but, as with our Sr°° example, this 
does not indicate that the ocean is the source of 
sulfur. 


Sincerely yours, 
CHRISTIAN E. JUNGE 


Air Force Cambridge Research Center 
L. G. Hanscom Field, Bedford, Mass. 
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Reply 


Dear Sir, 

The main question raised by Dr. Junge is no 
doubt of fundamental importance for the correct 
interpretation of chemistry data on precipitation, 
though his analysis perhaps can be extended some- 
what. Firstly, as most precipitation comes from dis- 
turbances which are moving comparatively slowly 
we can regard precipitation as a quasi-stationary 
state where horizontal convergence of moisture is 
steadily replacing the moisture precipitated. Neglect- 
ing contribution to the precipitation by sub-cloud 
air is probably permissible. Provided also that most 
of the precipitation is due to accretion of cloud drop- 
lets into larger aggregates, i.e. condensation on ice 
nuclei is quantitatively unimportant, the equation 
given by Dr. Junge is, of course, correct. But K in 
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his formula is an instantaneous value so it has to be 
averaged and the way it is done by Dr. Junge may 
be too simple. 

Denoting by m the amount of a substance precipi- 
tated, by R the rate of precipitation, and using in- 
stead of ı/l a specific volume v (volume per gram 
liquid water) we can write 


16 
m= / eRvcdt (2) 


fo) 


from which also the average weighted concentra- 
tion can be derived. The integration is, of course, 
carried out only during precipitation. 

Now we can assume ¢ to be constant, a reasonable 
assumption for most compounds of interest. In or- 
der to carry out the integration we use a common 
method writing R=R+R’, v=v+v" and c TA 
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Now an average of the sum of products like R’a’ 
can be written roro, + T where r is the correlation 
coefficient between R and v and or and 6, are stand- 
ard deviations of R and v respectively. Then we 
get for the integral 
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This looks complicated but, including the parenthe- 


sis into the constant we get 


= 


m=s R-v-ec-T (4) 


Now &’ depends on the correlations between and 
the relative variations of R, v and c. It may be con- 
stant within a smaller geographical unit but we do 
not know for certain. It may be different for mari- 
time compared to continental environments. We 
may assume it to be constant, however, but it may 
be less than ¢ or larger than ¢ depending on whether 
correlations are negative or positive. 

From eq. (4) we then get the average weighted 
concentration ¢, in precipitation 


D=E VC (5) 


and if v is constant in space ¢, will be entirely depend- 

‘ent on ¢. So when interpreting his maps of concen- 
trations in precipitation Dr. Junge assumes both v and 
&’ to be essentially constant over the whole area. Or, 
better, their variation does not obscure the linear re- 
lation between c, and ¢ too much. 

If R should be inversely related to v which means 
that the precipitation rate is proportional to the lig- 
uid water content in clouds, R :v should be essen- 
tially constant in a given area. As most precipitation 
comes from large scale disturbances, T over a longer 
time would also be constant. In addition, if ¢’ can 
be regarded as constant m should be proportional to 
€. This is then, in effect, the assumption I have made 
in my paper when interpreting total deposition data. 
If this assumption is correct, Dr. Junge’s assumption 
must be wrong and vice versa. 

Now, the rate of precipitation in large scale dis- 
turbances depends on the degree of horizontal con- 
vergence which can be amplified through mountains 
and ground friction. It depends also on the convec- 
tive processes going on in the clouds but, as these 
depend upon the release of latent heat through con- 
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densation they must depend on the liquid water con- 
tent. A lower value of v thus increases the convective 
activity and it is well known that this also causes an 
entrainment of upper, drier air into the clouds. The 
average concentration of sea salt is therefore more 
or less inversely proportional to the depth of the 
clouds and the same will be true more or less for 
any element originating at the earth’s surface which 
does not spend too long time in the atmosphere. 
Even for elements with a relatively long residence 
time, like sulphur, it will still hold close to the source 
areas. Hence, for elements originating at the earth’s 
surface one can expect R to be negatively correlated 
to c and possibly also to v. Then m will be propor- 
tional not to ¢ but to the total amount per unit area 
of the element. The concentration in precipitation 
will, however, be inversely related to the rate of 
precipitation. 


If, however, a chemical element should originate 


at the tropopause as in the case of Sr one would ex- 
pect something entirely different, namely that R and 
c should be positively correlated, especially as we 


know that the concentration of Sr increases with 
altitude. The effect of this positive correlation may 
be partly offset by a negative correlation between R 
and v. In this case v -¢ tends to be constant and m 
well largely depend on the amount of precipitation. 
This is an interesting point and it is rather significant 
that at least in Scandinavia the pattern of precipitat- 
ed radioactive debris does not follow that of precip- 
itated excess sulphur. In fact, nitrate-nitrogen is the 
only compound that bears any kind of resemblance 


90 . L'un { 
to Sr in Scandinavia, suggesting, at least partly, a 
stratospheric origin. 


It is thus seen that the interpretation of chemical 
data on precipitation is a rather complex problem 
and, until more direct evidence is available it is a 
matter of intuition and, perhaps, taste in which way 
they are interpreted. 


Dr. Junge also believes that the excess Na and K 
computed from this data for places close to the 
coasts are unreliable because they are obtained as dif- 
ferences between large numbers. However, as they 
are of the same sign, they cannot be due to random 
errors in the analyses. There are consequently only 
two alternatives left, one is that they are real and 
the other that they are due to systematic errors in 
the analyses. 


Dr. Junge also states that there is no basis for the 
assumption that the ionic composition of sea spray 
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particles deviates much from that of sea water. I 
have discussed some possibilities in my paper and I 
have the opinion that there is no basis for excluding 
the assumption that the ionic composition of sea 
spray particles deviates much from that of sea water. 

Finally Dr. Junge brings forth two references 
which he considers should have been included, and 
he is correct in that respect. Kalle, whom he refers 
to first used the data from the first Scandinavian 
network published by EMANUELSSON, ERIKSSON & 
EGNER in 1954. He too computed the excess Na, K 
and Ca in the precipitation at the different stations 
as I have done and plotted them as percentage of the 
total excess on the familiar triangular diagram to- 
gether with the average composition of various 
rocks with respect to these three components. 

He finds that the excess Na, K and Ca in precipi- 
tation falls outside the region of composition of 
various rocks. Now, this can be expected as any 
terrestrial matter must originate at the surface of 
soils and soils differ significantly in their composition 
from rocks. Kalle is also aware of this and makes 
some corrections on the excess to bring it in line 
with the composition of rocks. This is, of course, 
useless as rocks in general contribute nothing to the 
atmosphere, only the soils which are covering the 
rocks. Now, it seems also evident that the more sol- 
uble parts of the soil, at least in the areas which 
contribute to the atmosphere are dominated by the 
chemical composition of precipitation but not in a 
simple way. 

It is interesting though distressing to note Kalle’s 
opinion on the question of the origin of the excess 
components in precipitation. He states: «Vieles 
scheint dafür zu sprechen, dass der bisher unerklär- 
liche Vorgang der Jonentrennung bei der Nieder- 
schlagsbildung in der Atmosphäre durch den Ver- 


mischungsvorgang ozeanischer und litosphärischer 


«Some Statistical Aspects of the Dynamical Processes 
of Growth and Occlusion in Simple Baroclinic Models”, 
by Philip D. Thompson, Rossby Memorial Volume, 
pp. 350—358, The Rockefeller Institute Press, 1959. 


Dear Sir: 


In a recent paper by P. D. THOMPSON (1959) 
appearing in the Rossby Memorial Volume, a 
statistical dynamical approach was used to investigate 
the “growth” and “occlusion” in the 2-level 
baroclinic model. The main conclusions from the 
investigation are that 
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Salzanteile über den Kontinenten vorgetäuscht ist 
und dass es nicht nötig ist, zu völlig neuartigen, zum 
Teil reichlich hypothetischen Erklärungsversuchen 
seine Zuflucht nehmen zu müssen.” Principally, it is 
fair to disbelieve new theories and hypothesis, it is 
even necessary in science but one should never state 
categorically without documented arguments that 
there is no need for such theories and hypothesis. 

As for Meetham’s work I do not share Dr. Junge’s 
belief that industry is the dominant source of sulphur 
over England. I have also in my paper given some 
thought to this question. Both Dr. Junge and I agree 
that the residence time of sulphur in the atmosphere 
is appreciable, perhaps 2—3 weeks. Now, everybody 
acquainted to the British climate will have difficul- 
ties comprehending how the English industrially re- 
leased sulphur can be kept within the boundaries of 
England for that long. 

Meetham also computes a residence time of the 
industrial sulphur with respect to absorption by the 
ground and finds that it is somewhere between 5 and 
16 hours. If it were that short then one could dis- 
regard industrial sulphur as contributing to the aver- 
age concentration of sulphur in the atmosphere. 
However, apart from these points my view is that 
Meetham’s work is a significant contribution but his 
interpretation of the data may need some revision. 


Sincerely yours, 


ERIK ERIKSSON 


International Meteorological Institute 
Lindhagensgatan 124 
Stockholm K 


REPEREN CE: 


EMANUELSSON, A., ERIKSSON, E., and EGNÉR, H., 1954: 
Composition of atmospheric precipitation in Sweden, 
Tellus, 6, pp. 261—267. 


1. ‘ The simple 2-level model has a selective pre- 
ference for disturbances in which the temperature 
pattern lags behind the pressure pattern; its dynam- 
ical properties are such that disturbances in which 
the temperature field precedes the pressure field 
die out and remain undetected, while disturbances 
in which the temperature field lags behind are 
amplified.” 


‘ 


2. “... through the very same mechanism by 
which ‘out-of-phase’ disturbances grow, the tem- 
perature field is gradually brought into phase with 
the pressure field.” 
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Some reservation is apparently taken to the last 
conclusion in a note at the end of the paper, in 
which it is pointed out that there exists a stable 
stationary phase lag between the temperature and 
pressure fields provided the wave number exceeds 
a certain critical value. 

The conclusions quoted above from Thompson’s 
paper are based on a theory in which the velocity 
field is independent of the north-south direction 
and, in addition, that the amplitudes of the disturb- 
ances are small. Due to these assumptions one 
would expect to find an agreement between Thomp- 
son’s conclusions and the complete solution of the 
equations for the 2-level model, which can be 
solved in this particular case as demonstrated by 
Ocura (1957). 

From Ogura’s solutions it is first of all seen that 
there is no “growth” and no “occlusion” in the 
stable case. The amplitudes will have characteristic 
periods of fluctuation depending on the vertical 
windshear and the wave-number as seen from 
Table 1 in Ogura’s paper. The phase-angle will 
vary monotonically in the stable case as seen from 
his figure 2. 

In the unstable case Ogura has only illustrated 
one case, in which the disturbance at the upper 
level initially lags a quarter of a wave-length be- 
hind the disturbance at the lower level. It is found 
that the amplitude increases with time at both 
levels (fig. 3), and that the phase-angle approaches 
a limiting phase-angle depending on the vertical 
windshear and the wave-number (fig. 2). A closer 
inspection of the solution shows that the phase- 
angle independent of the initial state quite rapidly 
will approach the limiting phase-angle, which 
always is such that the disturbance at the upper 
level is lagging behind the disturbance at the lower 
level, or in Thompson’s formulation: the tempera- 
ture field is lagging behind the pressure field. 
With regard to the amplitudes of the disturbances 
it is found that they may decrease initially if the 
temperature field precedes the pressure fields, but 
that amplification sets in when the phase-angle has 
changed in such a way that the temperature field 
is lagging behind the pressure field. 

The 2-level model has therefore no selective 
preference for disturbances in which the temperature 
pattern lags behind the pressure pattern in the sense 
that disturbances with the opposite phase lag will 
die out. It is rather a result of the development 
of the flow pattern in the unstable case that the 
disturbances will have the temperature field lagging 
behind the pressure field. Another indication that 
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the 2-level model has no selective preference is the 
fact that the instability criterion for this model 
only depends on the wave-number, the vertical 
windshear, the static stability and the Coriolis pa- 
rameter. It has thus no reference to the initial ar- 
rangement of the temperature field relative to the 
pressure field. 

Concluding the comparison between Ogura’s 
and Thompson’s conclusions it is evident that the 
development of the flow in predictions with the 
2-level model, at least in the case where the flow is 
independent of the north-south direction, does not 
approach a state of quasi-barotropy, and that Thomp- 
son’s equations for “growth” and “occlusion” at 
best can give only the initial changes of the amplitude 
and phase-angle of the disturbance. Thompson’s 
conclusions are further in disagreement with recent 
calculations of the conversion between potential 
and kinetic energy based on the 2-level model as 
performed by SALTZMAN (1959) and the present 
writer (1959), in which it is shown that the conver- 
sion at any time is from potential to kinetic energy. 
This result can only be explained if the temperature 
field systematically lags behind the pressure field. 

Very truly yours, 
A. WrN-NIELSEN 


Joint Numerical Weather Prediction Unit, 
Suitland, Maryland 


Reply 
Dearssin: 


Several points of criticism implied (but charitably 
unstated) in Dr. Wiin-Nielsen’s letter I find to be 
very well taken. I have been uneasy about these 
points for some time and welcome this opportunity 
to set them straight. 

The conclusions quoted by Dr. Wiin-Nielsen, 
taken literally and as they stand, are correct only 
to the extent that they correctly relate the rate of 
growth to the current phase-lag; the error lay in 
misinterpreting my own equations and, in par- 
ticular, in supposing that the stationary states of 
amplitude and phase-lag are asymptotic states, in 
the sense that the system inevitably tends toward 
those states. In fact, I could not have resolved the 
latter question, simply because my formulation of 
the problem was then incomplete. 

A mathematically complete description of the 
behavior of the linear two-level model in question 
(sinusoidal perturbations of small amplitude, super- 


Tellus XII (1960), 4 


DEL TERS? TO DHE EDITOR 


posed on a horizontally uniform westerly motion, 
with velocity independent of latitude) is contained 
in the following equations: 


dR in aA 
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in which R is the ratio of the amplitudes at the 
upper and lower levels, respectively; À is the lag 
between the pressure perturbations at the two 
levels (positive for westward tilt); & is about the 
wave-number for minimum stability; U, and U, 
are the equilibrium values of the zonal wind at the 
lower and upper levels, respectively; and ß is 
the Rossby parameter. These equations state how 
the changes of amplitude and phase-lag are related 
to the amplitude and phase-lag themselves. 

To illustrate the nature of my earlier error, let us 
consider states in the neighborhood of the state 
(R =1, A=o). According to Eqs. (1) and (2), 
the latter state is a stationary one, but it is not an 
asymptotic state, simply because Eqs. (1) and (2) for 
small departures (R’ and 4’) from that stationary 
state reduce to: 


d2R’ pare 
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Thus, if the system lies near the state (R=1, À =o), 
it will oscillate around that state with frequency B/«, 
passing back and forth between positive and nega- 
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tive phase-lags. Dr. Wiin-Nielsen is quite correct 
in stating that the system will not invariably tend 
toward a quasi-barotropic state; it may, however, 
oscillate around a state of zero phase-lag, without 
regard to the vertical wind shear. Otherwise, a more 
complete investigation of Eqs. (1) and (2) leads 
to results that are in accord with Ogura’s conclusions 
for a case of large phase-lag. 

Some months before receiving Dr. Wiin-Nielsen’s 
letter, my uneasiness about these questions led me to 
undertake a study of the behavior of a linear 
baroclinic model with continuous vertical structure. 
In summary, I have derived a simultaneous system 
of two equations, analogous to Eqs. (1) and (2), 
in which the only unknowns are the vertical tilt 
and amplitude of the pressure perturbation (both 
considered as functions of time and elevation). 
Analysis of those equations indicates that a baroclinic 
model with continuous vertical structure does possess 
an asymptotic state—i.e., a state that is stationary 
in a special sense and which is stable in that the 
system always tends toward that state. Briefly 
described, this state is one in which the pressure 
perturbation has a definite vertical tilt toward the 
west, and in which amplification continues at all 
levels throughout the troposphere. These results 
will be reported in a forth-coming issue of Tellus. 

In short, far from emphasizing the virtues of the 
2-level model, I should stress its shortcomings. 
Although it may predict growth-rates that are 
initially correct, its mechanism of “occlusion” does 
not correctly simulate the atmosphere’s. For this 
reason alone, we must now study baroclinic models 
in which there is more freedom of interaction 
between the motions at different elevations. 


Sincerely yours, 
Pre D. THOMPSON 


International Meteorological Institute 
Lindhagensgatan 124 
Stockholm K 
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Erratum 


In Volume 12, Number 3 in the article “The daily variation of the cosmic ray nucleonic component 
at Murchison Bay and Uppsala” by A. E. SanpstrOM, E. Dyrıng and S. LINDGREN, the following correc- 
tion should be made: page 346, column 2, under Acknowledgements, line 8: change National to Natural. 
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